Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "NonlinearSystemBase.h"
11 : #include "AuxiliarySystem.h"
12 : #include "Problem.h"
13 : #include "FEProblem.h"
14 : #include "MooseVariableFE.h"
15 : #include "MooseVariableScalar.h"
16 : #include "PetscSupport.h"
17 : #include "Factory.h"
18 : #include "ParallelUniqueId.h"
19 : #include "ThreadedElementLoop.h"
20 : #include "MaterialData.h"
21 : #include "ComputeResidualThread.h"
22 : #include "ComputeResidualAndJacobianThread.h"
23 : #include "ComputeFVFluxThread.h"
24 : #include "ComputeJacobianThread.h"
25 : #include "ComputeJacobianForScalingThread.h"
26 : #include "ComputeFullJacobianThread.h"
27 : #include "ComputeJacobianBlocksThread.h"
28 : #include "ComputeDiracThread.h"
29 : #include "ComputeElemDampingThread.h"
30 : #include "ComputeNodalDampingThread.h"
31 : #include "ComputeNodalKernelsThread.h"
32 : #include "ComputeNodalKernelBcsThread.h"
33 : #include "ComputeNodalKernelJacobiansThread.h"
34 : #include "ComputeNodalKernelBCJacobiansThread.h"
35 : #include "TimeKernel.h"
36 : #include "BoundaryCondition.h"
37 : #include "DirichletBCBase.h"
38 : #include "NodalBCBase.h"
39 : #include "IntegratedBCBase.h"
40 : #include "DGKernel.h"
41 : #include "InterfaceKernelBase.h"
42 : #include "ElementDamper.h"
43 : #include "NodalDamper.h"
44 : #include "GeneralDamper.h"
45 : #include "DisplacedProblem.h"
46 : #include "NearestNodeLocator.h"
47 : #include "PenetrationLocator.h"
48 : #include "NodalConstraint.h"
49 : #include "NodeFaceConstraint.h"
50 : #include "NodeElemConstraintBase.h"
51 : #include "MortarConstraint.h"
52 : #include "ElemElemConstraint.h"
53 : #include "ScalarKernelBase.h"
54 : #include "Parser.h"
55 : #include "Split.h"
56 : #include "FieldSplitPreconditioner.h"
57 : #include "MooseMesh.h"
58 : #include "MooseUtils.h"
59 : #include "MooseApp.h"
60 : #include "NodalKernelBase.h"
61 : #include "DiracKernelBase.h"
62 : #include "TimeIntegrator.h"
63 : #include "Predictor.h"
64 : #include "Assembly.h"
65 : #include "ElementPairLocator.h"
66 : #include "ODETimeKernel.h"
67 : #include "AllLocalDofIndicesThread.h"
68 : #include "FloatingPointExceptionGuard.h"
69 : #include "ADKernel.h"
70 : #include "ADDirichletBCBase.h"
71 : #include "Moose.h"
72 : #include "ConsoleStream.h"
73 : #include "MooseError.h"
74 : #include "FVElementalKernel.h"
75 : #include "FVScalarLagrangeMultiplierConstraint.h"
76 : #include "FVBoundaryScalarLagrangeMultiplierConstraint.h"
77 : #include "FVFluxKernel.h"
78 : #include "FVScalarLagrangeMultiplierInterface.h"
79 : #include "UserObject.h"
80 : #include "OffDiagonalScalingMatrix.h"
81 : #include "HDGKernel.h"
82 :
83 : // libMesh
84 : #include "libmesh/nonlinear_solver.h"
85 : #include "libmesh/quadrature_gauss.h"
86 : #include "libmesh/dense_vector.h"
87 : #include "libmesh/boundary_info.h"
88 : #include "libmesh/petsc_matrix.h"
89 : #include "libmesh/petsc_vector.h"
90 : #include "libmesh/petsc_nonlinear_solver.h"
91 : #include "libmesh/numeric_vector.h"
92 : #include "libmesh/mesh.h"
93 : #include "libmesh/dense_subvector.h"
94 : #include "libmesh/dense_submatrix.h"
95 : #include "libmesh/dof_map.h"
96 : #include "libmesh/sparse_matrix.h"
97 : #include "libmesh/petsc_matrix.h"
98 : #include "libmesh/default_coupling.h"
99 : #include "libmesh/diagonal_matrix.h"
100 : #include "libmesh/fe_interface.h"
101 : #include "libmesh/petsc_solver_exception.h"
102 :
103 : #include <ios>
104 :
105 : #include "petscsnes.h"
106 : #include <PetscDMMoose.h>
107 : EXTERN_C_BEGIN
108 : extern PetscErrorCode DMCreate_Moose(DM);
109 : EXTERN_C_END
110 :
111 : using namespace libMesh;
112 :
113 62186 : NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem,
114 : System & sys,
115 62186 : const std::string & name)
116 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
117 62186 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
118 62186 : _sys(sys),
119 62186 : _last_nl_rnorm(0.),
120 62186 : _current_nl_its(0),
121 62186 : _residual_ghosted(NULL),
122 62186 : _Re_time_tag(-1),
123 62186 : _Re_time(NULL),
124 62186 : _Re_non_time_tag(-1),
125 62186 : _Re_non_time(NULL),
126 62186 : _scalar_kernels(/*threaded=*/false),
127 62186 : _nodal_bcs(/*threaded=*/false),
128 62186 : _preset_nodal_bcs(/*threaded=*/false),
129 62186 : _ad_preset_nodal_bcs(/*threaded=*/false),
130 62186 : _splits(/*threaded=*/false),
131 62186 : _increment_vec(NULL),
132 62186 : _use_finite_differenced_preconditioner(false),
133 62186 : _fdcoloring(nullptr),
134 62186 : _fsp(nullptr),
135 62186 : _add_implicit_geometric_coupling_entries_to_jacobian(false),
136 62186 : _assemble_constraints_separately(false),
137 62186 : _need_residual_ghosted(false),
138 62186 : _debugging_residuals(false),
139 62186 : _doing_dg(false),
140 62186 : _n_iters(0),
141 62186 : _n_linear_iters(0),
142 62186 : _n_residual_evaluations(0),
143 62186 : _final_residual(0.),
144 62186 : _computing_pre_smo_residual(false),
145 62186 : _pre_smo_residual(0),
146 62186 : _initial_residual(0),
147 62186 : _use_pre_smo_residual(false),
148 62186 : _print_all_var_norms(false),
149 62186 : _has_save_in(false),
150 62186 : _has_diag_save_in(false),
151 62186 : _has_nodalbc_save_in(false),
152 62186 : _has_nodalbc_diag_save_in(false),
153 62186 : _computed_scaling(false),
154 62186 : _compute_scaling_once(true),
155 62186 : _resid_vs_jac_scaling_param(0),
156 62186 : _off_diagonals_in_auto_scaling(false),
157 435302 : _auto_scaling_initd(false)
158 : {
159 62186 : getResidualNonTimeVector();
160 : // Don't need to add the matrix - it already exists (for now)
161 62186 : _Ke_system_tag = _fe_problem.addMatrixTag("SYSTEM");
162 :
163 : // The time matrix tag is not normally used - but must be added to the system
164 : // in case it is so that objects can have 'time' in their matrix tags by default
165 62186 : _fe_problem.addMatrixTag("TIME");
166 :
167 62186 : _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
168 :
169 62186 : _sys.identify_variable_groups(_fe_problem.identifyVariableGroupsInNL());
170 :
171 62186 : if (!_fe_problem.defaultGhosting())
172 : {
173 62104 : auto & dof_map = _sys.get_dof_map();
174 62104 : dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting());
175 62104 : dof_map.set_implicit_neighbor_dofs(false);
176 : }
177 62186 : }
178 :
179 57816 : NonlinearSystemBase::~NonlinearSystemBase() = default;
180 :
181 : void
182 60284 : NonlinearSystemBase::preInit()
183 : {
184 60284 : SolverSystem::preInit();
185 :
186 60284 : if (_fe_problem.hasDampers())
187 159 : setupDampers();
188 :
189 60284 : if (_residual_copy.get())
190 0 : _residual_copy->init(_sys.n_dofs(), false, SERIAL);
191 60284 : }
192 :
193 : void
194 7573 : NonlinearSystemBase::reinitMortarFunctors()
195 : {
196 : // reinit is called on meshChanged() in FEProblemBase. We could implement meshChanged() instead.
197 : // Subdomains might have changed
198 7649 : for (auto & functor : _displaced_mortar_functors)
199 76 : functor.second.setupMortarMaterials();
200 7689 : for (auto & functor : _undisplaced_mortar_functors)
201 116 : functor.second.setupMortarMaterials();
202 7573 : }
203 :
204 : void
205 76 : NonlinearSystemBase::turnOffJacobian()
206 : {
207 76 : system().set_basic_system_only();
208 76 : nonlinearSolver()->jacobian = NULL;
209 76 : }
210 :
211 : void
212 59380 : NonlinearSystemBase::initialSetup()
213 : {
214 296900 : TIME_SECTION("nlInitialSetup", 2, "Setting Up Nonlinear System");
215 :
216 59380 : SolverSystem::initialSetup();
217 :
218 : {
219 296900 : TIME_SECTION("kernelsInitialSetup", 2, "Setting Up Kernels/BCs/Constraints");
220 :
221 124329 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
222 : {
223 64949 : _kernels.initialSetup(tid);
224 64949 : _nodal_kernels.initialSetup(tid);
225 64949 : _dirac_kernels.initialSetup(tid);
226 64949 : if (_doing_dg)
227 1200 : _dg_kernels.initialSetup(tid);
228 64949 : _interface_kernels.initialSetup(tid);
229 :
230 64949 : _element_dampers.initialSetup(tid);
231 64949 : _nodal_dampers.initialSetup(tid);
232 64949 : _integrated_bcs.initialSetup(tid);
233 :
234 64949 : if (_fe_problem.haveFV())
235 : {
236 4324 : std::vector<FVElementalKernel *> fv_elemental_kernels;
237 4324 : _fe_problem.theWarehouse()
238 8648 : .query()
239 4324 : .template condition<AttribSystem>("FVElementalKernel")
240 4324 : .template condition<AttribThread>(tid)
241 4324 : .queryInto(fv_elemental_kernels);
242 :
243 8312 : for (auto * fv_kernel : fv_elemental_kernels)
244 3988 : fv_kernel->initialSetup();
245 :
246 4324 : std::vector<FVFluxKernel *> fv_flux_kernels;
247 4324 : _fe_problem.theWarehouse()
248 8648 : .query()
249 4324 : .template condition<AttribSystem>("FVFluxKernel")
250 4324 : .template condition<AttribThread>(tid)
251 4324 : .queryInto(fv_flux_kernels);
252 :
253 9767 : for (auto * fv_kernel : fv_flux_kernels)
254 5443 : fv_kernel->initialSetup();
255 4324 : }
256 : }
257 :
258 59380 : _scalar_kernels.initialSetup();
259 59380 : _constraints.initialSetup();
260 59380 : _general_dampers.initialSetup();
261 59380 : _nodal_bcs.initialSetup();
262 59380 : }
263 :
264 : {
265 296900 : TIME_SECTION("mortarSetup", 2, "Initializing Mortar Interfaces");
266 :
267 118760 : auto create_mortar_functors = [this](const bool displaced)
268 : {
269 : // go over mortar interfaces and construct functors
270 118760 : const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
271 119780 : for (const auto & mortar_interface : mortar_interfaces)
272 : {
273 1020 : const auto primary_secondary_boundary_pair = mortar_interface.first;
274 1020 : if (!_constraints.hasActiveMortarConstraints(primary_secondary_boundary_pair, displaced))
275 39 : continue;
276 :
277 981 : const auto & mortar_generation_object = mortar_interface.second;
278 :
279 : auto & mortar_constraints =
280 981 : _constraints.getActiveMortarConstraints(primary_secondary_boundary_pair, displaced);
281 :
282 : auto & subproblem = displaced
283 222 : ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
284 1092 : : static_cast<SubProblem &>(_fe_problem);
285 :
286 981 : auto & mortar_functors =
287 981 : displaced ? _displaced_mortar_functors : _undisplaced_mortar_functors;
288 :
289 981 : mortar_functors.emplace(primary_secondary_boundary_pair,
290 1962 : ComputeMortarFunctor(mortar_constraints,
291 : mortar_generation_object,
292 : subproblem,
293 : _fe_problem,
294 : displaced,
295 981 : subproblem.assembly(0, number())));
296 : }
297 118760 : };
298 :
299 59380 : create_mortar_functors(false);
300 59380 : create_mortar_functors(true);
301 59380 : }
302 :
303 59380 : if (_automatic_scaling)
304 : {
305 445 : if (_off_diagonals_in_auto_scaling)
306 69 : _scaling_matrix = std::make_unique<OffDiagonalScalingMatrix<Number>>(_communicator);
307 : else
308 376 : _scaling_matrix = std::make_unique<DiagonalMatrix<Number>>(_communicator);
309 : }
310 :
311 59380 : if (_preconditioner)
312 14278 : _preconditioner->initialSetup();
313 59380 : }
314 :
315 : void
316 287469 : NonlinearSystemBase::timestepSetup()
317 : {
318 287469 : SolverSystem::timestepSetup();
319 :
320 600590 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
321 : {
322 313121 : _kernels.timestepSetup(tid);
323 313121 : _nodal_kernels.timestepSetup(tid);
324 313121 : _dirac_kernels.timestepSetup(tid);
325 313121 : if (_doing_dg)
326 1775 : _dg_kernels.timestepSetup(tid);
327 313121 : _interface_kernels.timestepSetup(tid);
328 313121 : _element_dampers.timestepSetup(tid);
329 313121 : _nodal_dampers.timestepSetup(tid);
330 313121 : _integrated_bcs.timestepSetup(tid);
331 :
332 313121 : if (_fe_problem.haveFV())
333 : {
334 17435 : std::vector<FVFluxBC *> bcs;
335 17435 : _fe_problem.theWarehouse()
336 34870 : .query()
337 17435 : .template condition<AttribSystem>("FVFluxBC")
338 17435 : .template condition<AttribThread>(tid)
339 17435 : .queryInto(bcs);
340 :
341 17435 : std::vector<FVInterfaceKernel *> iks;
342 17435 : _fe_problem.theWarehouse()
343 34870 : .query()
344 17435 : .template condition<AttribSystem>("FVInterfaceKernel")
345 17435 : .template condition<AttribThread>(tid)
346 17435 : .queryInto(iks);
347 :
348 17435 : std::vector<FVFluxKernel *> kernels;
349 17435 : _fe_problem.theWarehouse()
350 34870 : .query()
351 17435 : .template condition<AttribSystem>("FVFluxKernel")
352 17435 : .template condition<AttribThread>(tid)
353 17435 : .queryInto(kernels);
354 :
355 34439 : for (auto * bc : bcs)
356 17004 : bc->timestepSetup();
357 17833 : for (auto * ik : iks)
358 398 : ik->timestepSetup();
359 38598 : for (auto * kernel : kernels)
360 21163 : kernel->timestepSetup();
361 17435 : }
362 : }
363 287469 : _scalar_kernels.timestepSetup();
364 287469 : _constraints.timestepSetup();
365 287469 : _general_dampers.timestepSetup();
366 287469 : _nodal_bcs.timestepSetup();
367 287469 : }
368 :
369 : void
370 1923833 : NonlinearSystemBase::customSetup(const ExecFlagType & exec_type)
371 : {
372 1923833 : SolverSystem::customSetup(exec_type);
373 :
374 4020605 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
375 : {
376 2096772 : _kernels.customSetup(exec_type, tid);
377 2096772 : _nodal_kernels.customSetup(exec_type, tid);
378 2096772 : _dirac_kernels.customSetup(exec_type, tid);
379 2096772 : if (_doing_dg)
380 12149 : _dg_kernels.customSetup(exec_type, tid);
381 2096772 : _interface_kernels.customSetup(exec_type, tid);
382 2096772 : _element_dampers.customSetup(exec_type, tid);
383 2096772 : _nodal_dampers.customSetup(exec_type, tid);
384 2096772 : _integrated_bcs.customSetup(exec_type, tid);
385 :
386 2096772 : if (_fe_problem.haveFV())
387 : {
388 128830 : std::vector<FVFluxBC *> bcs;
389 128830 : _fe_problem.theWarehouse()
390 257660 : .query()
391 128830 : .template condition<AttribSystem>("FVFluxBC")
392 128830 : .template condition<AttribThread>(tid)
393 128830 : .queryInto(bcs);
394 :
395 128830 : std::vector<FVInterfaceKernel *> iks;
396 128830 : _fe_problem.theWarehouse()
397 257660 : .query()
398 128830 : .template condition<AttribSystem>("FVInterfaceKernel")
399 128830 : .template condition<AttribThread>(tid)
400 128830 : .queryInto(iks);
401 :
402 128830 : std::vector<FVFluxKernel *> kernels;
403 128830 : _fe_problem.theWarehouse()
404 257660 : .query()
405 128830 : .template condition<AttribSystem>("FVFluxKernel")
406 128830 : .template condition<AttribThread>(tid)
407 128830 : .queryInto(kernels);
408 :
409 232205 : for (auto * bc : bcs)
410 103375 : bc->customSetup(exec_type);
411 173049 : for (auto * ik : iks)
412 44219 : ik->customSetup(exec_type);
413 306443 : for (auto * kernel : kernels)
414 177613 : kernel->customSetup(exec_type);
415 128830 : }
416 : }
417 1923833 : _scalar_kernels.customSetup(exec_type);
418 1923833 : _constraints.customSetup(exec_type);
419 1923833 : _general_dampers.customSetup(exec_type);
420 1923833 : _nodal_bcs.customSetup(exec_type);
421 1923833 : }
422 :
423 : void
424 343727 : NonlinearSystemBase::setupDM()
425 : {
426 343727 : if (_fsp)
427 129 : _fsp->setupDM();
428 343727 : }
429 :
430 : void
431 82636 : NonlinearSystemBase::addKernel(const std::string & kernel_name,
432 : const std::string & name,
433 : InputParameters & parameters)
434 : {
435 172768 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
436 : {
437 : // Create the kernel object via the factory and add to warehouse
438 : std::shared_ptr<KernelBase> kernel =
439 90332 : _factory.create<KernelBase>(kernel_name, name, parameters, tid);
440 90156 : _kernels.addObject(kernel, tid);
441 90136 : postAddResidualObject(*kernel);
442 : // Add to theWarehouse, a centralized storage for all moose objects
443 90132 : _fe_problem.theWarehouse().add(kernel);
444 90132 : }
445 :
446 82436 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
447 130 : _has_save_in = true;
448 82436 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
449 104 : _has_diag_save_in = true;
450 82436 : }
451 :
452 : void
453 808 : NonlinearSystemBase::addHDGKernel(const std::string & kernel_name,
454 : const std::string & name,
455 : InputParameters & parameters)
456 : {
457 1621 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
458 : {
459 : // Create the kernel object via the factory and add to warehouse
460 813 : auto kernel = _factory.create<HDGKernel>(kernel_name, name, parameters, tid);
461 813 : _kernels.addObject(kernel, tid);
462 813 : _hybridized_kernels.addObject(kernel, tid);
463 : // Add to theWarehouse, a centralized storage for all moose objects
464 813 : _fe_problem.theWarehouse().add(kernel);
465 813 : postAddResidualObject(*kernel);
466 813 : }
467 808 : }
468 :
469 : void
470 570 : NonlinearSystemBase::addNodalKernel(const std::string & kernel_name,
471 : const std::string & name,
472 : InputParameters & parameters)
473 : {
474 1241 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
475 : {
476 : // Create the kernel object via the factory and add to the warehouse
477 : std::shared_ptr<NodalKernelBase> kernel =
478 671 : _factory.create<NodalKernelBase>(kernel_name, name, parameters, tid);
479 671 : _nodal_kernels.addObject(kernel, tid);
480 : // Add to theWarehouse, a centralized storage for all moose objects
481 671 : _fe_problem.theWarehouse().add(kernel);
482 671 : postAddResidualObject(*kernel);
483 671 : }
484 :
485 570 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
486 0 : _has_save_in = true;
487 570 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
488 0 : _has_diag_save_in = true;
489 570 : }
490 :
491 : void
492 1551 : NonlinearSystemBase::addScalarKernel(const std::string & kernel_name,
493 : const std::string & name,
494 : InputParameters & parameters)
495 : {
496 : std::shared_ptr<ScalarKernelBase> kernel =
497 1551 : _factory.create<ScalarKernelBase>(kernel_name, name, parameters);
498 1543 : postAddResidualObject(*kernel);
499 : // Add to theWarehouse, a centralized storage for all moose objects
500 1543 : _fe_problem.theWarehouse().add(kernel);
501 1543 : _scalar_kernels.addObject(kernel);
502 1543 : }
503 :
504 : void
505 82315 : NonlinearSystemBase::addBoundaryCondition(const std::string & bc_name,
506 : const std::string & name,
507 : InputParameters & parameters)
508 : {
509 : // ThreadID
510 82315 : THREAD_ID tid = 0;
511 :
512 : // Create the object
513 : std::shared_ptr<BoundaryCondition> bc =
514 82315 : _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
515 82262 : postAddResidualObject(*bc);
516 :
517 : // Active BoundaryIDs for the object
518 82262 : const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
519 82262 : auto bc_var = dynamic_cast<const MooseVariableFieldBase *>(&bc->variable());
520 82262 : _vars[tid].addBoundaryVar(boundary_ids, bc_var);
521 :
522 : // Cast to the various types of BCs
523 82262 : std::shared_ptr<NodalBCBase> nbc = std::dynamic_pointer_cast<NodalBCBase>(bc);
524 82262 : std::shared_ptr<IntegratedBCBase> ibc = std::dynamic_pointer_cast<IntegratedBCBase>(bc);
525 :
526 : // NodalBCBase
527 82262 : if (nbc)
528 : {
529 73216 : if (nbc->checkNodalVar() && !nbc->variable().isNodal())
530 4 : mooseError("Trying to use nodal boundary condition '",
531 4 : nbc->name(),
532 : "' on a non-nodal variable '",
533 4 : nbc->variable().name(),
534 : "'.");
535 :
536 73212 : _nodal_bcs.addObject(nbc);
537 : // Add to theWarehouse, a centralized storage for all moose objects
538 73212 : _fe_problem.theWarehouse().add(nbc);
539 73212 : _vars[tid].addBoundaryVars(boundary_ids, nbc->getCoupledVars());
540 :
541 73212 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
542 99 : _has_nodalbc_save_in = true;
543 73212 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
544 21 : _has_nodalbc_diag_save_in = true;
545 :
546 : // DirichletBCs that are preset
547 73212 : std::shared_ptr<DirichletBCBase> dbc = std::dynamic_pointer_cast<DirichletBCBase>(bc);
548 73212 : if (dbc && dbc->preset())
549 63039 : _preset_nodal_bcs.addObject(dbc);
550 :
551 73212 : std::shared_ptr<ADDirichletBCBase> addbc = std::dynamic_pointer_cast<ADDirichletBCBase>(bc);
552 73212 : if (addbc && addbc->preset())
553 1636 : _ad_preset_nodal_bcs.addObject(addbc);
554 73212 : }
555 :
556 : // IntegratedBCBase
557 9046 : else if (ibc)
558 : {
559 9046 : _integrated_bcs.addObject(ibc, tid);
560 : // Add to theWarehouse, a centralized storage for all moose objects
561 9046 : _fe_problem.theWarehouse().add(ibc);
562 9046 : _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
563 :
564 9046 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
565 86 : _has_save_in = true;
566 9046 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
567 60 : _has_diag_save_in = true;
568 :
569 9722 : for (tid = 1; tid < libMesh::n_threads(); tid++)
570 : {
571 : // Create the object
572 676 : bc = _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
573 :
574 : // Give users opportunity to set some parameters
575 676 : postAddResidualObject(*bc);
576 :
577 : // Active BoundaryIDs for the object
578 676 : const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
579 676 : _vars[tid].addBoundaryVar(boundary_ids, bc_var);
580 :
581 676 : ibc = std::static_pointer_cast<IntegratedBCBase>(bc);
582 :
583 676 : _integrated_bcs.addObject(ibc, tid);
584 676 : _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
585 : }
586 : }
587 :
588 : else
589 0 : mooseError("Unknown BoundaryCondition type for object named ", bc->name());
590 82258 : }
591 :
592 : void
593 1727 : NonlinearSystemBase::addConstraint(const std::string & c_name,
594 : const std::string & name,
595 : InputParameters & parameters)
596 : {
597 1727 : std::shared_ptr<Constraint> constraint = _factory.create<Constraint>(c_name, name, parameters);
598 1719 : _constraints.addObject(constraint);
599 1719 : postAddResidualObject(*constraint);
600 :
601 1719 : if (constraint && constraint->addCouplingEntriesToJacobian())
602 1706 : addImplicitGeometricCouplingEntriesToJacobian(true);
603 1719 : }
604 :
605 : void
606 915 : NonlinearSystemBase::addDiracKernel(const std::string & kernel_name,
607 : const std::string & name,
608 : InputParameters & parameters)
609 : {
610 1906 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
611 : {
612 : std::shared_ptr<DiracKernelBase> kernel =
613 999 : _factory.create<DiracKernelBase>(kernel_name, name, parameters, tid);
614 991 : postAddResidualObject(*kernel);
615 991 : _dirac_kernels.addObject(kernel, tid);
616 : // Add to theWarehouse, a centralized storage for all moose objects
617 991 : _fe_problem.theWarehouse().add(kernel);
618 991 : }
619 907 : }
620 :
621 : void
622 1403 : NonlinearSystemBase::addDGKernel(std::string dg_kernel_name,
623 : const std::string & name,
624 : InputParameters & parameters)
625 : {
626 2926 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
627 : {
628 1523 : auto dg_kernel = _factory.create<DGKernelBase>(dg_kernel_name, name, parameters, tid);
629 1523 : _dg_kernels.addObject(dg_kernel, tid);
630 : // Add to theWarehouse, a centralized storage for all moose objects
631 1523 : _fe_problem.theWarehouse().add(dg_kernel);
632 1523 : postAddResidualObject(*dg_kernel);
633 1523 : }
634 :
635 1403 : _doing_dg = true;
636 :
637 1403 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
638 52 : _has_save_in = true;
639 1403 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
640 39 : _has_diag_save_in = true;
641 1403 : }
642 :
643 : void
644 907 : NonlinearSystemBase::addInterfaceKernel(std::string interface_kernel_name,
645 : const std::string & name,
646 : InputParameters & parameters)
647 : {
648 1889 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
649 : {
650 : std::shared_ptr<InterfaceKernelBase> interface_kernel =
651 982 : _factory.create<InterfaceKernelBase>(interface_kernel_name, name, parameters, tid);
652 982 : postAddResidualObject(*interface_kernel);
653 :
654 982 : const std::set<BoundaryID> & boundary_ids = interface_kernel->boundaryIDs();
655 982 : auto ik_var = dynamic_cast<const MooseVariableFieldBase *>(&interface_kernel->variable());
656 982 : _vars[tid].addBoundaryVar(boundary_ids, ik_var);
657 :
658 982 : _interface_kernels.addObject(interface_kernel, tid);
659 : // Add to theWarehouse, a centralized storage for all moose objects
660 982 : _fe_problem.theWarehouse().add(interface_kernel);
661 982 : _vars[tid].addBoundaryVars(boundary_ids, interface_kernel->getCoupledVars());
662 982 : }
663 907 : }
664 :
665 : void
666 177 : NonlinearSystemBase::addDamper(const std::string & damper_name,
667 : const std::string & name,
668 : InputParameters & parameters)
669 : {
670 328 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
671 : {
672 218 : std::shared_ptr<Damper> damper = _factory.create<Damper>(damper_name, name, parameters, tid);
673 :
674 : // Attempt to cast to the damper types
675 218 : std::shared_ptr<ElementDamper> ed = std::dynamic_pointer_cast<ElementDamper>(damper);
676 218 : std::shared_ptr<NodalDamper> nd = std::dynamic_pointer_cast<NodalDamper>(damper);
677 218 : std::shared_ptr<GeneralDamper> gd = std::dynamic_pointer_cast<GeneralDamper>(damper);
678 :
679 218 : if (gd)
680 : {
681 67 : _general_dampers.addObject(gd);
682 67 : break; // not threaded
683 : }
684 151 : else if (ed)
685 92 : _element_dampers.addObject(ed, tid);
686 59 : else if (nd)
687 59 : _nodal_dampers.addObject(nd, tid);
688 : else
689 0 : mooseError("Invalid damper type");
690 419 : }
691 177 : }
692 :
693 : void
694 485 : NonlinearSystemBase::addSplit(const std::string & split_name,
695 : const std::string & name,
696 : InputParameters & parameters)
697 : {
698 485 : std::shared_ptr<Split> split = _factory.create<Split>(split_name, name, parameters);
699 485 : _splits.addObject(split);
700 : // Add to theWarehouse, a centralized storage for all moose objects
701 485 : _fe_problem.theWarehouse().add(split);
702 485 : }
703 :
704 : std::shared_ptr<Split>
705 485 : NonlinearSystemBase::getSplit(const std::string & name)
706 : {
707 485 : return _splits.getActiveObject(name);
708 : }
709 :
710 : bool
711 312713 : NonlinearSystemBase::shouldEvaluatePreSMOResidual() const
712 : {
713 312713 : if (_fe_problem.solverParams(number())._type == Moose::ST_LINEAR)
714 10638 : return false;
715 :
716 : // The legacy behavior (#10464) _always_ performs the pre-SMO residual evaluation
717 : // regardless of whether it is needed.
718 : //
719 : // This is not ideal and has been fixed by #23472. This legacy option ensures a smooth transition
720 : // to the new behavior. Modules and Apps that want to migrate to the new behavior should set this
721 : // parameter to false.
722 302075 : if (_app.parameters().get<bool>("use_legacy_initial_residual_evaluation_behavior"))
723 0 : return true;
724 :
725 302075 : return _use_pre_smo_residual;
726 : }
727 :
728 : Real
729 386861 : NonlinearSystemBase::referenceResidual() const
730 : {
731 386861 : return usePreSMOResidual() ? preSMOResidual() : initialResidual();
732 : }
733 :
734 : Real
735 171 : NonlinearSystemBase::preSMOResidual() const
736 : {
737 171 : if (!shouldEvaluatePreSMOResidual())
738 0 : mooseError("pre-SMO residual is requested but not evaluated.");
739 :
740 171 : return _pre_smo_residual;
741 : }
742 :
743 : Real
744 386973 : NonlinearSystemBase::initialResidual() const
745 : {
746 386973 : return _initial_residual;
747 : }
748 :
749 : void
750 310100 : NonlinearSystemBase::setInitialResidual(Real r)
751 : {
752 310100 : _initial_residual = r;
753 310100 : }
754 :
755 : void
756 0 : NonlinearSystemBase::zeroVectorForResidual(const std::string & vector_name)
757 : {
758 0 : for (unsigned int i = 0; i < _vecs_to_zero_for_residual.size(); ++i)
759 0 : if (vector_name == _vecs_to_zero_for_residual[i])
760 0 : return;
761 :
762 0 : _vecs_to_zero_for_residual.push_back(vector_name);
763 : }
764 :
765 : void
766 131 : NonlinearSystemBase::computeResidualTag(NumericVector<Number> & residual, TagID tag_id)
767 : {
768 131 : _nl_vector_tags.clear();
769 131 : _nl_vector_tags.insert(tag_id);
770 131 : _nl_vector_tags.insert(residualVectorTag());
771 :
772 131 : associateVectorToTag(residual, residualVectorTag());
773 :
774 131 : computeResidualTags(_nl_vector_tags);
775 :
776 131 : disassociateVectorFromTag(residual, residualVectorTag());
777 131 : }
778 :
779 : void
780 0 : NonlinearSystemBase::computeResidual(NumericVector<Number> & residual, TagID tag_id)
781 : {
782 0 : mooseDeprecated(" Please use computeResidualTag");
783 :
784 0 : computeResidualTag(residual, tag_id);
785 0 : }
786 :
787 : void
788 3293480 : NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
789 : {
790 : parallel_object_only();
791 :
792 9880440 : TIME_SECTION("nl::computeResidualTags", 5);
793 :
794 3293480 : _fe_problem.setCurrentNonlinearSystem(number());
795 3293480 : _fe_problem.setCurrentlyComputingResidual(true);
796 :
797 3293480 : bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
798 :
799 3293480 : _n_residual_evaluations++;
800 :
801 : // not suppose to do anythin on matrix
802 3293480 : deactivateAllMatrixTags();
803 :
804 3293480 : FloatingPointExceptionGuard fpe_guard(_app);
805 :
806 3293480 : for (const auto & numeric_vec : _vecs_to_zero_for_residual)
807 0 : if (hasVector(numeric_vec))
808 : {
809 0 : NumericVector<Number> & vec = getVector(numeric_vec);
810 0 : vec.close();
811 0 : vec.zero();
812 : }
813 :
814 : try
815 : {
816 3293480 : zeroTaggedVectors(tags);
817 3293480 : computeResidualInternal(tags);
818 3293315 : closeTaggedVectors(tags);
819 :
820 3293315 : if (required_residual)
821 : {
822 3261433 : auto & residual = getVector(residualVectorTag());
823 3261433 : if (!_time_integrators.empty())
824 : {
825 5738717 : for (auto & ti : _time_integrators)
826 2877784 : ti->postResidual(residual);
827 : }
828 : else
829 400500 : residual += *_Re_non_time;
830 3261433 : residual.close();
831 : }
832 3293315 : if (_fe_problem.computingScalingResidual())
833 : // We don't want to do nodal bcs or anything else
834 50 : return;
835 :
836 3293265 : computeNodalBCs(tags);
837 3293265 : closeTaggedVectors(tags);
838 :
839 : // If we are debugging residuals we need one more assignment to have the ghosted copy up to
840 : // date
841 3293265 : if (_need_residual_ghosted && _debugging_residuals && required_residual)
842 : {
843 2010 : auto & residual = getVector(residualVectorTag());
844 :
845 2010 : *_residual_ghosted = residual;
846 2010 : _residual_ghosted->close();
847 : }
848 : // Need to close and update the aux system in case residuals were saved to it.
849 3293265 : if (_has_nodalbc_save_in)
850 171 : _fe_problem.getAuxiliarySystem().solution().close();
851 3293265 : if (hasSaveIn())
852 309 : _fe_problem.getAuxiliarySystem().update();
853 : }
854 121 : catch (MooseException & e)
855 : {
856 : // The buck stops here, we have already handled the exception by
857 : // calling stopSolve(), it is now up to PETSc to return a
858 : // "diverged" reason during the next solve.
859 121 : }
860 :
861 : // not supposed to do anything on matrix
862 3293386 : activateAllMatrixTags();
863 :
864 3293386 : _fe_problem.setCurrentlyComputingResidual(false);
865 3293486 : }
866 :
867 : void
868 3467 : NonlinearSystemBase::computeResidualAndJacobianTags(const std::set<TagID> & vector_tags,
869 : const std::set<TagID> & matrix_tags)
870 : {
871 : const bool required_residual =
872 3467 : vector_tags.find(residualVectorTag()) == vector_tags.end() ? false : true;
873 :
874 : try
875 : {
876 3467 : zeroTaggedVectors(vector_tags);
877 3467 : computeResidualAndJacobianInternal(vector_tags, matrix_tags);
878 3467 : closeTaggedVectors(vector_tags);
879 3467 : closeTaggedMatrices(matrix_tags);
880 :
881 3467 : if (required_residual)
882 : {
883 3467 : auto & residual = getVector(residualVectorTag());
884 3467 : if (!_time_integrators.empty())
885 : {
886 3818 : for (auto & ti : _time_integrators)
887 1909 : ti->postResidual(residual);
888 : }
889 : else
890 1558 : residual += *_Re_non_time;
891 3467 : residual.close();
892 : }
893 :
894 3467 : computeNodalBCsResidualAndJacobian();
895 3467 : closeTaggedVectors(vector_tags);
896 3467 : closeTaggedMatrices(matrix_tags);
897 : }
898 0 : catch (MooseException & e)
899 : {
900 : // The buck stops here, we have already handled the exception by
901 : // calling stopSolve(), it is now up to PETSc to return a
902 : // "diverged" reason during the next solve.
903 0 : }
904 3467 : }
905 :
906 : void
907 258647 : NonlinearSystemBase::onTimestepBegin()
908 : {
909 518012 : for (auto & ti : _time_integrators)
910 259365 : ti->preSolve();
911 258647 : if (_predictor.get())
912 387 : _predictor->timestepSetup();
913 258647 : }
914 :
915 : void
916 313848 : NonlinearSystemBase::setInitialSolution()
917 : {
918 313848 : deactivateAllMatrixTags();
919 :
920 313848 : NumericVector<Number> & initial_solution(solution());
921 313848 : if (_predictor.get())
922 : {
923 387 : if (_predictor->shouldApply())
924 : {
925 1115 : TIME_SECTION("applyPredictor", 2, "Applying Predictor");
926 :
927 223 : _predictor->apply(initial_solution);
928 223 : _fe_problem.predictorCleanup(initial_solution);
929 223 : }
930 : else
931 164 : _console << " Skipping predictor this step" << std::endl;
932 : }
933 :
934 : // do nodal BC
935 : {
936 1569240 : TIME_SECTION("initialBCs", 2, "Applying BCs To Initial Condition");
937 :
938 313848 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
939 15392653 : for (const auto & bnode : bnd_nodes)
940 : {
941 15078805 : BoundaryID boundary_id = bnode->_bnd_id;
942 15078805 : Node * node = bnode->_node;
943 :
944 15078805 : if (node->processor_id() == processor_id())
945 : {
946 : // reinit variables in nodes
947 11559322 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
948 :
949 11559322 : if (_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
950 : {
951 3634329 : const auto & preset_bcs = _preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
952 7923110 : for (const auto & preset_bc : preset_bcs)
953 4288781 : preset_bc->computeValue(initial_solution);
954 : }
955 11559322 : if (_ad_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
956 : {
957 36720 : const auto & preset_bcs_res = _ad_preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
958 77278 : for (const auto & preset_bc : preset_bcs_res)
959 40558 : preset_bc->computeValue(initial_solution);
960 : }
961 : }
962 : }
963 313848 : }
964 :
965 313848 : _sys.solution->close();
966 313848 : update();
967 :
968 : // Set constraint secondary values
969 313848 : setConstraintSecondaryValues(initial_solution, false);
970 :
971 313848 : if (_fe_problem.getDisplacedProblem())
972 35761 : setConstraintSecondaryValues(initial_solution, true);
973 313848 : }
974 :
975 : void
976 121 : NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
977 : {
978 121 : _predictor = predictor;
979 121 : }
980 :
981 : void
982 6188103 : NonlinearSystemBase::subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
983 : {
984 6188103 : SolverSystem::subdomainSetup();
985 :
986 6188103 : _kernels.subdomainSetup(subdomain, tid);
987 6188103 : _nodal_kernels.subdomainSetup(subdomain, tid);
988 6188103 : _element_dampers.subdomainSetup(subdomain, tid);
989 6188103 : _nodal_dampers.subdomainSetup(subdomain, tid);
990 6188103 : }
991 :
992 : NumericVector<Number> &
993 31305 : NonlinearSystemBase::getResidualTimeVector()
994 : {
995 31305 : if (!_Re_time)
996 : {
997 31199 : _Re_time_tag = _fe_problem.addVectorTag("TIME");
998 :
999 : // Most applications don't need the expense of ghosting
1000 31199 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1001 31199 : _Re_time = &addVector(_Re_time_tag, false, ptype);
1002 : }
1003 106 : else if (_need_residual_ghosted && _Re_time->type() == PARALLEL)
1004 : {
1005 0 : const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
1006 :
1007 : // If an application changes its mind, the libMesh API lets us
1008 : // change the vector.
1009 0 : _Re_time = &system().add_vector(vector_name, false, GHOSTED);
1010 0 : }
1011 :
1012 31305 : return *_Re_time;
1013 : }
1014 :
1015 : NumericVector<Number> &
1016 93698 : NonlinearSystemBase::getResidualNonTimeVector()
1017 : {
1018 93698 : if (!_Re_non_time)
1019 : {
1020 62186 : _Re_non_time_tag = _fe_problem.addVectorTag("NONTIME");
1021 :
1022 : // Most applications don't need the expense of ghosting
1023 62186 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1024 62186 : _Re_non_time = &addVector(_Re_non_time_tag, false, ptype);
1025 : }
1026 31512 : else if (_need_residual_ghosted && _Re_non_time->type() == PARALLEL)
1027 : {
1028 0 : const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
1029 :
1030 : // If an application changes its mind, the libMesh API lets us
1031 : // change the vector.
1032 0 : _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
1033 0 : }
1034 :
1035 93698 : return *_Re_non_time;
1036 : }
1037 :
1038 : NumericVector<Number> &
1039 0 : NonlinearSystemBase::residualVector(TagID tag)
1040 : {
1041 0 : mooseDeprecated("Please use getVector()");
1042 0 : switch (tag)
1043 : {
1044 0 : case 0:
1045 0 : return getResidualNonTimeVector();
1046 :
1047 0 : case 1:
1048 0 : return getResidualTimeVector();
1049 :
1050 0 : default:
1051 0 : mooseError("The required residual vector is not available");
1052 : }
1053 : }
1054 :
1055 : void
1056 23409 : NonlinearSystemBase::enforceNodalConstraintsResidual(NumericVector<Number> & residual)
1057 : {
1058 23409 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1059 23409 : residual.close();
1060 23409 : if (_constraints.hasActiveNodalConstraints())
1061 : {
1062 602 : const auto & ncs = _constraints.getActiveNodalConstraints();
1063 1204 : for (const auto & nc : ncs)
1064 : {
1065 602 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1066 602 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1067 :
1068 602 : if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1069 : {
1070 602 : _fe_problem.reinitNodes(primary_node_ids, tid);
1071 602 : _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1072 602 : nc->computeResidual(residual);
1073 : }
1074 : }
1075 602 : _fe_problem.addCachedResidualDirectly(residual, tid);
1076 602 : residual.close();
1077 : }
1078 23409 : }
1079 :
1080 : void
1081 3667 : NonlinearSystemBase::enforceNodalConstraintsJacobian()
1082 : {
1083 3667 : if (!hasMatrix(systemMatrixTag()))
1084 0 : mooseError(" A system matrix is required");
1085 :
1086 3667 : auto & jacobian = getMatrix(systemMatrixTag());
1087 3667 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1088 :
1089 3667 : if (_constraints.hasActiveNodalConstraints())
1090 : {
1091 121 : const auto & ncs = _constraints.getActiveNodalConstraints();
1092 242 : for (const auto & nc : ncs)
1093 : {
1094 121 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1095 121 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1096 :
1097 121 : if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1098 : {
1099 121 : _fe_problem.reinitNodes(primary_node_ids, tid);
1100 121 : _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1101 121 : nc->computeJacobian(jacobian);
1102 : }
1103 : }
1104 121 : _fe_problem.addCachedJacobian(tid);
1105 : }
1106 3667 : }
1107 :
1108 : void
1109 108178 : NonlinearSystemBase::reinitNodeFace(const Node & secondary_node,
1110 : const BoundaryID secondary_boundary,
1111 : const PenetrationInfo & info,
1112 : const bool displaced)
1113 : {
1114 115408 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1115 165882 : : static_cast<SubProblem &>(_fe_problem);
1116 :
1117 108178 : const Elem * primary_elem = info._elem;
1118 108178 : unsigned int primary_side = info._side_num;
1119 108178 : std::vector<Point> points;
1120 108178 : points.push_back(info._closest_point);
1121 :
1122 : // *These next steps MUST be done in this order!*
1123 : // ADL: This is a Chesterton's fence situation. I don't know which calls exactly the above comment
1124 : // is referring to. If I had to guess I would guess just the reinitNodeFace and prepareAssembly
1125 : // calls since the former will size the variable's dof indices and then the latter will resize the
1126 : // residual/Jacobian based off the variable's cached dof indices size
1127 :
1128 : // This reinits the variables that exist on the secondary node
1129 108178 : _fe_problem.reinitNodeFace(&secondary_node, secondary_boundary, 0);
1130 :
1131 : // This will set aside residual and jacobian space for the variables that have dofs on
1132 : // the secondary node
1133 108178 : _fe_problem.prepareAssembly(0);
1134 :
1135 108178 : _fe_problem.setNeighborSubdomainID(primary_elem, 0);
1136 :
1137 : //
1138 : // Reinit material on undisplaced mesh
1139 : //
1140 :
1141 : const Elem * const undisplaced_primary_elem =
1142 108178 : displaced ? _mesh.elemPtr(primary_elem->id()) : primary_elem;
1143 : const Point undisplaced_primary_physical_point =
1144 0 : [&points, displaced, primary_elem, undisplaced_primary_elem]()
1145 : {
1146 108178 : if (displaced)
1147 : {
1148 : const Point reference_point =
1149 57704 : FEMap::inverse_map(primary_elem->dim(), primary_elem, points[0]);
1150 57704 : return FEMap::map(primary_elem->dim(), undisplaced_primary_elem, reference_point);
1151 : }
1152 : else
1153 : // If our penetration locator is on the reference mesh, then our undisplaced
1154 : // physical point is simply the point coming from the penetration locator
1155 50474 : return points[0];
1156 108178 : }();
1157 :
1158 216356 : _fe_problem.reinitNeighborPhys(
1159 : undisplaced_primary_elem, primary_side, {undisplaced_primary_physical_point}, 0);
1160 : // Stateful material properties are only initialized for neighbor material data for internal faces
1161 : // for discontinuous Galerkin methods or for conforming interfaces for interface kernels. We don't
1162 : // have either of those use cases here where we likely have disconnected meshes
1163 108178 : _fe_problem.reinitMaterialsNeighbor(primary_elem->subdomain_id(), 0, /*swap_stateful=*/false);
1164 :
1165 : // Reinit points for constraint enforcement
1166 108178 : if (displaced)
1167 57704 : subproblem.reinitNeighborPhys(primary_elem, primary_side, points, 0);
1168 108178 : }
1169 :
1170 : void
1171 349609 : NonlinearSystemBase::setConstraintSecondaryValues(NumericVector<Number> & solution, bool displaced)
1172 : {
1173 :
1174 : if (displaced)
1175 : mooseAssert(_fe_problem.getDisplacedProblem(),
1176 : "If we're calling this method with displaced = true, then we better well have a "
1177 : "displaced problem");
1178 71522 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1179 385370 : : static_cast<SubProblem &>(_fe_problem);
1180 349609 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1181 :
1182 349609 : bool constraints_applied = false;
1183 :
1184 409830 : for (const auto & it : penetration_locators)
1185 : {
1186 60221 : PenetrationLocator & pen_loc = *(it.second);
1187 :
1188 60221 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1189 :
1190 60221 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1191 60221 : BoundaryID primary_boundary = pen_loc._primary_boundary;
1192 :
1193 60221 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1194 : {
1195 : const auto & constraints =
1196 2977 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1197 2977 : std::unordered_set<unsigned int> needed_mat_props;
1198 5954 : for (const auto & constraint : constraints)
1199 : {
1200 2977 : const auto & mp_deps = constraint->getMatPropDependencies();
1201 2977 : needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
1202 : }
1203 2977 : _fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/0);
1204 :
1205 42140 : for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1206 : {
1207 39163 : dof_id_type secondary_node_num = secondary_nodes[i];
1208 39163 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1209 :
1210 39163 : if (secondary_node.processor_id() == processor_id())
1211 : {
1212 30944 : if (pen_loc._penetration_info[secondary_node_num])
1213 : {
1214 30944 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1215 :
1216 30944 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1217 :
1218 61888 : for (const auto & nfc : constraints)
1219 : {
1220 30944 : if (nfc->isExplicitConstraint())
1221 28816 : continue;
1222 : // Return if this constraint does not correspond to the primary-secondary pair
1223 : // prepared by the outer loops.
1224 : // This continue statement is required when, e.g. one secondary surface constrains
1225 : // more than one primary surface.
1226 4256 : if (nfc->secondaryBoundary() != secondary_boundary ||
1227 2128 : nfc->primaryBoundary() != primary_boundary)
1228 0 : continue;
1229 :
1230 2128 : if (nfc->shouldApply())
1231 : {
1232 2128 : constraints_applied = true;
1233 2128 : nfc->computeSecondaryValue(solution);
1234 : }
1235 :
1236 2128 : if (nfc->hasWritableCoupledVariables())
1237 : {
1238 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1239 0 : for (auto * var : nfc->getWritableCoupledVariables())
1240 : {
1241 0 : if (var->isNodalDefined())
1242 0 : var->insert(_fe_problem.getAuxiliarySystem().solution());
1243 : }
1244 0 : }
1245 : }
1246 : }
1247 : }
1248 : }
1249 2977 : }
1250 : }
1251 :
1252 : // go over NodeELemConstraints
1253 349609 : std::set<dof_id_type> unique_secondary_node_ids;
1254 :
1255 785458 : for (const auto & secondary_id : _mesh.meshSubdomains())
1256 : {
1257 1165798 : for (const auto & primary_id : _mesh.meshSubdomains())
1258 : {
1259 729949 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1260 : {
1261 : const auto & constraints =
1262 180 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1263 :
1264 : // get unique set of ids of all nodes on current block
1265 180 : unique_secondary_node_ids.clear();
1266 180 : const MeshBase & meshhelper = _mesh.getMesh();
1267 360 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1268 17340 : meshhelper.active_subdomain_elements_end(secondary_id)))
1269 : {
1270 49140 : for (auto & n : elem->node_ref_range())
1271 40740 : unique_secondary_node_ids.insert(n.id());
1272 180 : }
1273 :
1274 14040 : for (auto secondary_node_id : unique_secondary_node_ids)
1275 : {
1276 13860 : Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1277 :
1278 : // check if secondary node is on current processor
1279 13860 : if (secondary_node.processor_id() == processor_id())
1280 : {
1281 : // This reinits the variables that exist on the secondary node
1282 11088 : _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
1283 :
1284 : // This will set aside residual and jacobian space for the variables that have dofs
1285 : // on the secondary node
1286 11088 : _fe_problem.prepareAssembly(0);
1287 :
1288 23584 : for (const auto & nec : constraints)
1289 : {
1290 12496 : if (nec->shouldApply())
1291 : {
1292 5680 : constraints_applied = true;
1293 5680 : nec->computeSecondaryValue(solution);
1294 : }
1295 : }
1296 : }
1297 : }
1298 : }
1299 : }
1300 : }
1301 :
1302 : // See if constraints were applied anywhere
1303 349609 : _communicator.max(constraints_applied);
1304 :
1305 349609 : if (constraints_applied)
1306 : {
1307 719 : solution.close();
1308 719 : update();
1309 : }
1310 349609 : }
1311 :
1312 : void
1313 28287 : NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool displaced)
1314 : {
1315 : // Make sure the residual is in a good state
1316 28287 : residual.close();
1317 :
1318 : if (displaced)
1319 : mooseAssert(_fe_problem.getDisplacedProblem(),
1320 : "If we're calling this method with displaced = true, then we better well have a "
1321 : "displaced problem");
1322 9756 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1323 33165 : : static_cast<SubProblem &>(_fe_problem);
1324 28287 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1325 :
1326 : bool constraints_applied;
1327 28287 : bool residual_has_inserted_values = false;
1328 28287 : if (!_assemble_constraints_separately)
1329 28287 : constraints_applied = false;
1330 52353 : for (const auto & it : penetration_locators)
1331 : {
1332 24066 : if (_assemble_constraints_separately)
1333 : {
1334 : // Reset the constraint_applied flag before each new constraint, as they need to be
1335 : // assembled separately
1336 0 : constraints_applied = false;
1337 : }
1338 24066 : PenetrationLocator & pen_loc = *(it.second);
1339 :
1340 24066 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1341 :
1342 24066 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1343 24066 : BoundaryID primary_boundary = pen_loc._primary_boundary;
1344 :
1345 24066 : bool has_writable_variables(false);
1346 :
1347 24066 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1348 : {
1349 : const auto & constraints =
1350 14004 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1351 :
1352 99926 : for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1353 : {
1354 85922 : dof_id_type secondary_node_num = secondary_nodes[i];
1355 85922 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1356 :
1357 85922 : if (secondary_node.processor_id() == processor_id())
1358 : {
1359 73050 : if (pen_loc._penetration_info[secondary_node_num])
1360 : {
1361 73050 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1362 :
1363 73050 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1364 :
1365 146100 : for (const auto & nfc : constraints)
1366 : {
1367 : // Return if this constraint does not correspond to the primary-secondary pair
1368 : // prepared by the outer loops.
1369 : // This continue statement is required when, e.g. one secondary surface constrains
1370 : // more than one primary surface.
1371 146100 : if (nfc->secondaryBoundary() != secondary_boundary ||
1372 73050 : nfc->primaryBoundary() != primary_boundary)
1373 0 : continue;
1374 :
1375 73050 : if (nfc->shouldApply())
1376 : {
1377 44234 : constraints_applied = true;
1378 44234 : nfc->computeResidual();
1379 :
1380 44234 : if (nfc->overwriteSecondaryResidual())
1381 : {
1382 : // The below will actually overwrite the residual for every single dof that
1383 : // lives on the node. We definitely don't want to do that!
1384 : // _fe_problem.setResidual(residual, 0);
1385 :
1386 44234 : const auto & secondary_var = nfc->variable();
1387 44234 : const auto & secondary_dofs = secondary_var.dofIndices();
1388 : mooseAssert(secondary_dofs.size() == secondary_var.count(),
1389 : "We are on a node so there should only be one dof per variable (for "
1390 : "an ArrayVariable we should have a number of dofs equal to the "
1391 : "number of components");
1392 :
1393 : // Assume that if the user is overwriting the secondary residual, then they are
1394 : // supplying residuals that do not correspond to their other physics
1395 : // (e.g. Kernels), hence we should not apply a scalingFactor that is normally
1396 : // based on the order of their other physics (e.g. Kernels)
1397 88468 : std::vector<Number> values = {nfc->secondaryResidual()};
1398 44234 : residual.insert(values, secondary_dofs);
1399 44234 : residual_has_inserted_values = true;
1400 44234 : }
1401 : else
1402 0 : _fe_problem.cacheResidual(0);
1403 44234 : _fe_problem.cacheResidualNeighbor(0);
1404 : }
1405 73050 : if (nfc->hasWritableCoupledVariables())
1406 : {
1407 28816 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1408 28816 : has_writable_variables = true;
1409 57632 : for (auto * var : nfc->getWritableCoupledVariables())
1410 : {
1411 28816 : if (var->isNodalDefined())
1412 28816 : var->insert(_fe_problem.getAuxiliarySystem().solution());
1413 : }
1414 28816 : }
1415 : }
1416 : }
1417 : }
1418 : }
1419 : }
1420 24066 : _communicator.max(has_writable_variables);
1421 :
1422 24066 : if (has_writable_variables)
1423 : {
1424 : // Explicit contact dynamic constraints write to auxiliary variables and update the old
1425 : // displacement solution on the constraint boundaries. Close solutions and update system
1426 : // accordingly.
1427 2401 : _fe_problem.getAuxiliarySystem().solution().close();
1428 2401 : _fe_problem.getAuxiliarySystem().system().update();
1429 2401 : solutionOld().close();
1430 : }
1431 :
1432 24066 : if (_assemble_constraints_separately)
1433 : {
1434 : // Make sure that secondary contribution to primary are assembled, and ghosts have been
1435 : // exchanged, as current primaries might become secondaries on next iteration and will need to
1436 : // contribute their former secondaries' contributions to the future primaries. See if
1437 : // constraints were applied anywhere
1438 0 : _communicator.max(constraints_applied);
1439 :
1440 0 : if (constraints_applied)
1441 : {
1442 : // If any of the above constraints inserted values in the residual, it needs to be
1443 : // assembled before adding the cached residuals below.
1444 0 : _communicator.max(residual_has_inserted_values);
1445 0 : if (residual_has_inserted_values)
1446 : {
1447 0 : residual.close();
1448 0 : residual_has_inserted_values = false;
1449 : }
1450 0 : _fe_problem.addCachedResidualDirectly(residual, 0);
1451 0 : residual.close();
1452 :
1453 0 : if (_need_residual_ghosted)
1454 0 : *_residual_ghosted = residual;
1455 : }
1456 : }
1457 : }
1458 28287 : if (!_assemble_constraints_separately)
1459 : {
1460 28287 : _communicator.max(constraints_applied);
1461 :
1462 28287 : if (constraints_applied)
1463 : {
1464 : // If any of the above constraints inserted values in the residual, it needs to be assembled
1465 : // before adding the cached residuals below.
1466 10864 : _communicator.max(residual_has_inserted_values);
1467 10864 : if (residual_has_inserted_values)
1468 10864 : residual.close();
1469 :
1470 10864 : _fe_problem.addCachedResidualDirectly(residual, 0);
1471 10864 : residual.close();
1472 :
1473 10864 : if (_need_residual_ghosted)
1474 10864 : *_residual_ghosted = residual;
1475 : }
1476 : }
1477 :
1478 : // go over element-element constraint interface
1479 28287 : THREAD_ID tid = 0;
1480 28287 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
1481 28287 : for (const auto & it : element_pair_locators)
1482 : {
1483 0 : ElementPairLocator & elem_pair_loc = *(it.second);
1484 :
1485 0 : if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
1486 : {
1487 : // ElemElemConstraint objects
1488 : const auto & _element_constraints =
1489 0 : _constraints.getActiveElemElemConstraints(it.first, displaced);
1490 :
1491 : // go over pair elements
1492 : const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
1493 0 : elem_pair_loc.getElemPairs();
1494 0 : for (const auto & pr : elem_pairs)
1495 : {
1496 0 : const Elem * elem1 = pr.first;
1497 0 : const Elem * elem2 = pr.second;
1498 :
1499 0 : if (elem1->processor_id() != processor_id())
1500 0 : continue;
1501 :
1502 0 : const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
1503 :
1504 : // for each element process constraints on the
1505 0 : for (const auto & ec : _element_constraints)
1506 : {
1507 0 : _fe_problem.setCurrentSubdomainID(elem1, tid);
1508 0 : subproblem.reinitElemPhys(elem1, info._elem1_constraint_q_point, tid);
1509 0 : _fe_problem.setNeighborSubdomainID(elem2, tid);
1510 0 : subproblem.reinitNeighborPhys(elem2, info._elem2_constraint_q_point, tid);
1511 :
1512 0 : ec->prepareShapes(ec->variable().number());
1513 0 : ec->prepareNeighborShapes(ec->variable().number());
1514 :
1515 0 : ec->reinit(info);
1516 0 : ec->computeResidual();
1517 0 : _fe_problem.cacheResidual(tid);
1518 0 : _fe_problem.cacheResidualNeighbor(tid);
1519 : }
1520 0 : _fe_problem.addCachedResidual(tid);
1521 : }
1522 : }
1523 : }
1524 :
1525 : // go over NodeElemConstraints
1526 28287 : std::set<dof_id_type> unique_secondary_node_ids;
1527 :
1528 28287 : constraints_applied = false;
1529 28287 : residual_has_inserted_values = false;
1530 28287 : bool has_writable_variables = false;
1531 97398 : for (const auto & secondary_id : _mesh.meshSubdomains())
1532 : {
1533 296138 : for (const auto & primary_id : _mesh.meshSubdomains())
1534 : {
1535 227027 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1536 : {
1537 : const auto & constraints =
1538 360 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1539 :
1540 : // get unique set of ids of all nodes on current block
1541 360 : unique_secondary_node_ids.clear();
1542 360 : const MeshBase & meshhelper = _mesh.getMesh();
1543 720 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1544 17880 : meshhelper.active_subdomain_elements_end(secondary_id)))
1545 : {
1546 98280 : for (auto & n : elem->node_ref_range())
1547 81480 : unique_secondary_node_ids.insert(n.id());
1548 360 : }
1549 :
1550 28080 : for (auto secondary_node_id : unique_secondary_node_ids)
1551 : {
1552 27720 : Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1553 : // check if secondary node is on current processor
1554 27720 : if (secondary_node.processor_id() == processor_id())
1555 : {
1556 : // This reinits the variables that exist on the secondary node
1557 22176 : _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
1558 :
1559 : // This will set aside residual and jacobian space for the variables that have dofs
1560 : // on the secondary node
1561 22176 : _fe_problem.prepareAssembly(0);
1562 :
1563 47168 : for (const auto & nec : constraints)
1564 : {
1565 24992 : if (nec->shouldApply())
1566 : {
1567 11360 : constraints_applied = true;
1568 11360 : nec->computeResidual();
1569 :
1570 11360 : if (nec->overwriteSecondaryResidual())
1571 : {
1572 0 : _fe_problem.setResidual(residual, 0);
1573 0 : residual_has_inserted_values = true;
1574 : }
1575 : else
1576 11360 : _fe_problem.cacheResidual(0);
1577 11360 : _fe_problem.cacheResidualNeighbor(0);
1578 : }
1579 24992 : if (nec->hasWritableCoupledVariables())
1580 : {
1581 352 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1582 352 : has_writable_variables = true;
1583 704 : for (auto * var : nec->getWritableCoupledVariables())
1584 : {
1585 352 : if (var->isNodalDefined())
1586 352 : var->insert(_fe_problem.getAuxiliarySystem().solution());
1587 : }
1588 352 : }
1589 : }
1590 22176 : _fe_problem.addCachedResidual(0);
1591 : }
1592 : }
1593 : }
1594 : }
1595 : }
1596 28287 : _communicator.max(constraints_applied);
1597 :
1598 28287 : if (constraints_applied)
1599 : {
1600 : // If any of the above constraints inserted values in the residual, it needs to be assembled
1601 : // before adding the cached residuals below.
1602 360 : _communicator.max(residual_has_inserted_values);
1603 360 : if (residual_has_inserted_values)
1604 0 : residual.close();
1605 :
1606 360 : _fe_problem.addCachedResidualDirectly(residual, 0);
1607 360 : residual.close();
1608 :
1609 360 : if (_need_residual_ghosted)
1610 360 : *_residual_ghosted = residual;
1611 : }
1612 28287 : _communicator.max(has_writable_variables);
1613 :
1614 28287 : if (has_writable_variables)
1615 : {
1616 : // Explicit contact dynamic constraints write to auxiliary variables and update the old
1617 : // displacement solution on the constraint boundaries. Close solutions and update system
1618 : // accordingly.
1619 20 : _fe_problem.getAuxiliarySystem().solution().close();
1620 20 : _fe_problem.getAuxiliarySystem().system().update();
1621 20 : solutionOld().close();
1622 : }
1623 :
1624 : // We may have additional tagged vectors that also need to be accumulated
1625 28287 : _fe_problem.addCachedResidual(0);
1626 28287 : }
1627 :
1628 : void
1629 4497 : NonlinearSystemBase::overwriteNodeFace(NumericVector<Number> & soln)
1630 : {
1631 : // Overwrite results from integrator in case we have explicit dynamics contact constraints
1632 4497 : auto & subproblem = _fe_problem.getDisplacedProblem()
1633 6898 : ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1634 6898 : : static_cast<SubProblem &>(_fe_problem);
1635 4497 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1636 :
1637 6898 : for (const auto & it : penetration_locators)
1638 : {
1639 2401 : PenetrationLocator & pen_loc = *(it.second);
1640 :
1641 2401 : const auto & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1642 2401 : const BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1643 2401 : const BoundaryID primary_boundary = pen_loc._primary_boundary;
1644 :
1645 2401 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, true))
1646 : {
1647 : const auto & constraints =
1648 2401 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, true);
1649 39217 : for (const auto i : index_range(secondary_nodes))
1650 : {
1651 36816 : const auto secondary_node_num = secondary_nodes[i];
1652 36816 : const Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1653 :
1654 36816 : if (secondary_node.processor_id() == processor_id())
1655 28816 : if (pen_loc._penetration_info[secondary_node_num])
1656 57632 : for (const auto & nfc : constraints)
1657 : {
1658 28816 : if (!nfc->isExplicitConstraint())
1659 0 : continue;
1660 :
1661 : // Return if this constraint does not correspond to the primary-secondary pair
1662 : // prepared by the outer loops.
1663 : // This continue statement is required when, e.g. one secondary surface constrains
1664 : // more than one primary surface.
1665 57632 : if (nfc->secondaryBoundary() != secondary_boundary ||
1666 28816 : nfc->primaryBoundary() != primary_boundary)
1667 0 : continue;
1668 :
1669 28816 : nfc->overwriteBoundaryVariables(soln, secondary_node);
1670 : }
1671 : }
1672 : }
1673 : }
1674 4497 : soln.close();
1675 4497 : }
1676 :
1677 : void
1678 3389487 : NonlinearSystemBase::residualSetup()
1679 : {
1680 10168461 : TIME_SECTION("residualSetup", 3);
1681 :
1682 3389487 : SolverSystem::residualSetup();
1683 :
1684 7082915 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1685 : {
1686 3693428 : _kernels.residualSetup(tid);
1687 3693428 : _nodal_kernels.residualSetup(tid);
1688 3693428 : _dirac_kernels.residualSetup(tid);
1689 3693428 : if (_doing_dg)
1690 35297 : _dg_kernels.residualSetup(tid);
1691 3693428 : _interface_kernels.residualSetup(tid);
1692 3693428 : _element_dampers.residualSetup(tid);
1693 3693428 : _nodal_dampers.residualSetup(tid);
1694 3693428 : _integrated_bcs.residualSetup(tid);
1695 : }
1696 3389487 : _scalar_kernels.residualSetup();
1697 3389487 : _constraints.residualSetup();
1698 3389487 : _general_dampers.residualSetup();
1699 3389487 : _nodal_bcs.residualSetup();
1700 :
1701 : // Avoid recursion
1702 3389487 : if (this == &_fe_problem.currentNonlinearSystem())
1703 3296947 : _fe_problem.residualSetup();
1704 3389487 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
1705 3389487 : }
1706 :
1707 : void
1708 3293480 : NonlinearSystemBase::computeResidualInternal(const std::set<TagID> & tags)
1709 : {
1710 : parallel_object_only();
1711 :
1712 9880440 : TIME_SECTION("computeResidualInternal", 3);
1713 :
1714 3293480 : residualSetup();
1715 :
1716 3293480 : const auto vector_tag_data = _fe_problem.getVectorTags(tags);
1717 :
1718 : // Residual contributions from UOs - for now this is used for ray tracing
1719 : // and ray kernels that contribute to the residual (think line sources)
1720 3293480 : std::vector<UserObject *> uos;
1721 3293480 : _fe_problem.theWarehouse()
1722 6586960 : .query()
1723 3293480 : .condition<AttribSystem>("UserObject")
1724 3293480 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1725 3293480 : .queryInto(uos);
1726 3293480 : for (auto & uo : uos)
1727 0 : uo->residualSetup();
1728 3293480 : for (auto & uo : uos)
1729 : {
1730 0 : uo->initialize();
1731 0 : uo->execute();
1732 0 : uo->finalize();
1733 : }
1734 :
1735 : // reinit scalar variables
1736 6882682 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1737 3589202 : _fe_problem.reinitScalars(tid);
1738 :
1739 : // residual contributions from the domain
1740 : PARALLEL_TRY
1741 : {
1742 9880440 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1743 :
1744 3293480 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
1745 :
1746 3293480 : ComputeResidualThread cr(_fe_problem, tags);
1747 3293480 : Threads::parallel_reduce(elem_range, cr);
1748 :
1749 : // We pass face information directly to FV residual objects for their evaluation. Consequently
1750 : // we must make sure to do separate threaded loops for 1) undisplaced face information objects
1751 : // and undisplaced residual objects and 2) displaced face information objects and displaced
1752 : // residual objects
1753 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
1754 3293456 : if (_fe_problem.haveFV())
1755 : {
1756 : ComputeFVFluxResidualThread<FVRange> fvr(
1757 150558 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
1758 150558 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
1759 150558 : Threads::parallel_reduce(faces, fvr);
1760 150558 : }
1761 3293448 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
1762 3293448 : displaced_problem && displaced_problem->haveFV())
1763 : {
1764 : ComputeFVFluxResidualThread<FVRange> fvr(
1765 0 : _fe_problem, this->number(), tags, /*on_displaced=*/true);
1766 0 : FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
1767 0 : displaced_problem->mesh().ownedFaceInfoEnd());
1768 0 : Threads::parallel_reduce(faces, fvr);
1769 3293448 : }
1770 :
1771 3293448 : unsigned int n_threads = libMesh::n_threads();
1772 6882610 : for (unsigned int i = 0; i < n_threads;
1773 : i++) // Add any cached residuals that might be hanging around
1774 3589162 : _fe_problem.addCachedResidual(i);
1775 3293448 : }
1776 3293448 : PARALLEL_CATCH;
1777 :
1778 : // residual contributions from the scalar kernels
1779 : PARALLEL_TRY
1780 : {
1781 : // do scalar kernels (not sure how to thread this)
1782 3293327 : if (_scalar_kernels.hasActiveObjects())
1783 : {
1784 221043 : TIME_SECTION("ScalarKernels", 3 /*, "Computing ScalarKernels"*/);
1785 :
1786 : MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
1787 : // This code should be refactored once we can do tags for scalar
1788 : // kernels
1789 : // Should redo this based on Warehouse
1790 73681 : if (!tags.size() || tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL))
1791 73660 : scalar_kernel_warehouse = &_scalar_kernels;
1792 21 : else if (tags.size() == 1)
1793 : scalar_kernel_warehouse =
1794 12 : &(_scalar_kernels.getVectorTagObjectWarehouse(*(tags.begin()), 0));
1795 : else
1796 : // scalar_kernels is not threading
1797 9 : scalar_kernel_warehouse = &(_scalar_kernels.getVectorTagsObjectWarehouse(tags, 0));
1798 :
1799 73681 : bool have_scalar_contributions = false;
1800 73681 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
1801 270946 : for (const auto & scalar_kernel : scalars)
1802 : {
1803 197265 : scalar_kernel->reinit();
1804 197265 : const std::vector<dof_id_type> & dof_indices = scalar_kernel->variable().dofIndices();
1805 197265 : const DofMap & dof_map = scalar_kernel->variable().dofMap();
1806 197265 : const dof_id_type first_dof = dof_map.first_dof();
1807 197265 : const dof_id_type end_dof = dof_map.end_dof();
1808 227734 : for (dof_id_type dof : dof_indices)
1809 : {
1810 197842 : if (dof >= first_dof && dof < end_dof)
1811 : {
1812 167373 : scalar_kernel->computeResidual();
1813 167373 : have_scalar_contributions = true;
1814 167373 : break;
1815 : }
1816 : }
1817 : }
1818 73681 : if (have_scalar_contributions)
1819 64576 : _fe_problem.addResidualScalar();
1820 73681 : }
1821 : }
1822 3293327 : PARALLEL_CATCH;
1823 :
1824 : // residual contributions from Block NodalKernels
1825 : PARALLEL_TRY
1826 : {
1827 3293327 : if (_nodal_kernels.hasActiveBlockObjects())
1828 : {
1829 46788 : TIME_SECTION("NodalKernels", 3 /*, "Computing NodalKernels"*/);
1830 :
1831 15596 : ComputeNodalKernelsThread cnk(_fe_problem, _nodal_kernels, tags);
1832 :
1833 15596 : const ConstNodeRange & range = _fe_problem.getCurrentAlgebraicNodeRange();
1834 :
1835 15596 : if (range.begin() != range.end())
1836 : {
1837 15596 : _fe_problem.reinitNode(*range.begin(), 0);
1838 :
1839 15596 : Threads::parallel_reduce(range, cnk);
1840 :
1841 15596 : unsigned int n_threads = libMesh::n_threads();
1842 38526 : for (unsigned int i = 0; i < n_threads;
1843 : i++) // Add any cached residuals that might be hanging around
1844 22930 : _fe_problem.addCachedResidual(i);
1845 : }
1846 15596 : }
1847 : }
1848 3293327 : PARALLEL_CATCH;
1849 :
1850 3293327 : if (_fe_problem.computingScalingResidual())
1851 : // We computed the volumetric objects. We can return now before we get into
1852 : // any strongly enforced constraint conditions or penalty-type objects
1853 : // (DGKernels, IntegratedBCs, InterfaceKernels, Constraints)
1854 50 : return;
1855 :
1856 : // residual contributions from boundary NodalKernels
1857 : PARALLEL_TRY
1858 : {
1859 3293277 : if (_nodal_kernels.hasActiveBoundaryObjects())
1860 : {
1861 3012 : TIME_SECTION("NodalKernelBCs", 3 /*, "Computing NodalKernelBCs"*/);
1862 :
1863 1004 : ComputeNodalKernelBcsThread cnk(_fe_problem, _nodal_kernels, tags);
1864 :
1865 1004 : const ConstBndNodeRange & bnd_node_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
1866 :
1867 1004 : Threads::parallel_reduce(bnd_node_range, cnk);
1868 :
1869 1004 : unsigned int n_threads = libMesh::n_threads();
1870 2094 : for (unsigned int i = 0; i < n_threads;
1871 : i++) // Add any cached residuals that might be hanging around
1872 1090 : _fe_problem.addCachedResidual(i);
1873 1004 : }
1874 : }
1875 3293277 : PARALLEL_CATCH;
1876 :
1877 3293277 : mortarConstraints(Moose::ComputeType::Residual, tags, {});
1878 :
1879 3293277 : if (_residual_copy.get())
1880 : {
1881 0 : _Re_non_time->close();
1882 0 : _Re_non_time->localize(*_residual_copy);
1883 : }
1884 :
1885 3293277 : if (_need_residual_ghosted)
1886 : {
1887 13234 : _Re_non_time->close();
1888 13234 : *_residual_ghosted = *_Re_non_time;
1889 13234 : _residual_ghosted->close();
1890 : }
1891 :
1892 3293277 : PARALLEL_TRY { computeDiracContributions(tags, false); }
1893 3293265 : PARALLEL_CATCH;
1894 :
1895 3293265 : if (_fe_problem._has_constraints)
1896 : {
1897 23409 : PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1898 23409 : PARALLEL_CATCH;
1899 23409 : _Re_non_time->close();
1900 : }
1901 :
1902 : // Add in Residual contributions from other Constraints
1903 3293265 : if (_fe_problem._has_constraints)
1904 : {
1905 : PARALLEL_TRY
1906 : {
1907 : // Undisplaced Constraints
1908 23409 : constraintResiduals(*_Re_non_time, false);
1909 :
1910 : // Displaced Constraints
1911 23409 : if (_fe_problem.getDisplacedProblem())
1912 4878 : constraintResiduals(*_Re_non_time, true);
1913 :
1914 23409 : if (_fe_problem.computingNonlinearResid())
1915 10705 : _constraints.residualEnd();
1916 : }
1917 23409 : PARALLEL_CATCH;
1918 23409 : _Re_non_time->close();
1919 : }
1920 :
1921 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
1922 : // counters
1923 3293265 : _app.solutionInvalidity().syncIteration();
1924 3293265 : _app.solutionInvalidity().solutionInvalidAccumulation();
1925 3293778 : }
1926 :
1927 : void
1928 3467 : NonlinearSystemBase::computeResidualAndJacobianInternal(const std::set<TagID> & vector_tags,
1929 : const std::set<TagID> & matrix_tags)
1930 : {
1931 10401 : TIME_SECTION("computeResidualAndJacobianInternal", 3);
1932 :
1933 : // Make matrix ready to use
1934 3467 : activateAllMatrixTags();
1935 :
1936 10401 : for (auto tag : matrix_tags)
1937 : {
1938 6934 : if (!hasMatrix(tag))
1939 3467 : continue;
1940 :
1941 3467 : auto & jacobian = getMatrix(tag);
1942 : // Necessary for speed
1943 3467 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
1944 : {
1945 3467 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
1946 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
1947 : PETSC_TRUE));
1948 3467 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
1949 0 : LibmeshPetscCall(
1950 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1951 3467 : if (_fe_problem.ignoreZerosInJacobian())
1952 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
1953 : MAT_IGNORE_ZERO_ENTRIES,
1954 : PETSC_TRUE));
1955 : }
1956 : }
1957 :
1958 3467 : residualSetup();
1959 :
1960 : // Residual contributions from UOs - for now this is used for ray tracing
1961 : // and ray kernels that contribute to the residual (think line sources)
1962 3467 : std::vector<UserObject *> uos;
1963 3467 : _fe_problem.theWarehouse()
1964 6934 : .query()
1965 3467 : .condition<AttribSystem>("UserObject")
1966 3467 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1967 3467 : .queryInto(uos);
1968 3467 : for (auto & uo : uos)
1969 0 : uo->residualSetup();
1970 3467 : for (auto & uo : uos)
1971 : {
1972 0 : uo->initialize();
1973 0 : uo->execute();
1974 0 : uo->finalize();
1975 : }
1976 :
1977 : // reinit scalar variables
1978 7273 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1979 3806 : _fe_problem.reinitScalars(tid);
1980 :
1981 : // residual contributions from the domain
1982 : PARALLEL_TRY
1983 : {
1984 10401 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1985 :
1986 3467 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
1987 :
1988 3467 : ComputeResidualAndJacobianThread crj(_fe_problem, vector_tags, matrix_tags);
1989 3467 : Threads::parallel_reduce(elem_range, crj);
1990 :
1991 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
1992 3467 : if (_fe_problem.haveFV())
1993 : {
1994 : ComputeFVFluxRJThread<FVRange> fvrj(
1995 1451 : _fe_problem, this->number(), vector_tags, matrix_tags, /*on_displaced=*/false);
1996 1451 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
1997 1451 : Threads::parallel_reduce(faces, fvrj);
1998 1451 : }
1999 3467 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
2000 3467 : displaced_problem && displaced_problem->haveFV())
2001 : {
2002 : ComputeFVFluxRJThread<FVRange> fvr(
2003 0 : _fe_problem, this->number(), vector_tags, matrix_tags, /*on_displaced=*/true);
2004 0 : FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
2005 0 : displaced_problem->mesh().ownedFaceInfoEnd());
2006 0 : Threads::parallel_reduce(faces, fvr);
2007 3467 : }
2008 :
2009 3467 : mortarConstraints(Moose::ComputeType::ResidualAndJacobian, vector_tags, matrix_tags);
2010 :
2011 3467 : unsigned int n_threads = libMesh::n_threads();
2012 7273 : for (unsigned int i = 0; i < n_threads;
2013 : i++) // Add any cached residuals that might be hanging around
2014 : {
2015 3806 : _fe_problem.addCachedResidual(i);
2016 3806 : _fe_problem.addCachedJacobian(i);
2017 : }
2018 3467 : }
2019 3467 : PARALLEL_CATCH;
2020 3467 : }
2021 :
2022 : void
2023 0 : NonlinearSystemBase::computeNodalBCs(NumericVector<Number> & residual)
2024 : {
2025 0 : _nl_vector_tags.clear();
2026 :
2027 0 : const auto & residual_vector_tags = _fe_problem.getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
2028 0 : for (const auto & residual_vector_tag : residual_vector_tags)
2029 0 : _nl_vector_tags.insert(residual_vector_tag._id);
2030 :
2031 0 : associateVectorToTag(residual, residualVectorTag());
2032 0 : computeNodalBCs(residual, _nl_vector_tags);
2033 0 : disassociateVectorFromTag(residual, residualVectorTag());
2034 0 : }
2035 :
2036 : void
2037 0 : NonlinearSystemBase::computeNodalBCs(NumericVector<Number> & residual, const std::set<TagID> & tags)
2038 : {
2039 0 : associateVectorToTag(residual, residualVectorTag());
2040 :
2041 0 : computeNodalBCs(tags);
2042 :
2043 0 : disassociateVectorFromTag(residual, residualVectorTag());
2044 0 : }
2045 :
2046 : void
2047 3293265 : NonlinearSystemBase::computeNodalBCs(const std::set<TagID> & tags)
2048 : {
2049 : // We need to close the diag_save_in variables on the aux system before NodalBCBases clear the
2050 : // dofs on boundary nodes
2051 3293265 : if (_has_save_in)
2052 309 : _fe_problem.getAuxiliarySystem().solution().close();
2053 :
2054 : // Select nodal kernels
2055 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2056 :
2057 3293265 : if (tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL) || !tags.size())
2058 3254434 : nbc_warehouse = &_nodal_bcs;
2059 38831 : else if (tags.size() == 1)
2060 18286 : nbc_warehouse = &(_nodal_bcs.getVectorTagObjectWarehouse(*(tags.begin()), 0));
2061 : else
2062 20545 : nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0));
2063 :
2064 : // Return early if there is no nodal kernel
2065 3293265 : if (!nbc_warehouse->size())
2066 420152 : return;
2067 :
2068 : PARALLEL_TRY
2069 : {
2070 2873113 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2071 :
2072 2873113 : if (!bnd_nodes.empty())
2073 : {
2074 8618406 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2075 :
2076 150801504 : for (const auto & bnode : bnd_nodes)
2077 : {
2078 147928702 : BoundaryID boundary_id = bnode->_bnd_id;
2079 147928702 : Node * node = bnode->_node;
2080 :
2081 261696152 : if (node->processor_id() == processor_id() &&
2082 113767450 : nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
2083 : {
2084 : // reinit variables in nodes
2085 58149792 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2086 :
2087 58149792 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2088 121908876 : for (const auto & nbc : bcs)
2089 63759084 : if (nbc->shouldApply())
2090 63757680 : nbc->computeResidual();
2091 : }
2092 : }
2093 2872802 : }
2094 : }
2095 2873113 : PARALLEL_CATCH;
2096 :
2097 2873113 : if (_Re_time)
2098 2551133 : _Re_time->close();
2099 2873113 : _Re_non_time->close();
2100 : }
2101 :
2102 : void
2103 3467 : NonlinearSystemBase::computeNodalBCsResidualAndJacobian()
2104 : {
2105 : PARALLEL_TRY
2106 : {
2107 3467 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2108 :
2109 3467 : if (!bnd_nodes.empty())
2110 : {
2111 10401 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2112 :
2113 86555 : for (const auto & bnode : bnd_nodes)
2114 : {
2115 83088 : BoundaryID boundary_id = bnode->_bnd_id;
2116 83088 : Node * node = bnode->_node;
2117 :
2118 83088 : if (node->processor_id() == processor_id())
2119 : {
2120 : // reinit variables in nodes
2121 62704 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2122 62704 : if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
2123 : {
2124 12501 : const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
2125 25002 : for (const auto & nbc : bcs)
2126 12501 : if (nbc->shouldApply())
2127 12501 : nbc->computeResidualAndJacobian();
2128 : }
2129 : }
2130 : }
2131 3467 : }
2132 : }
2133 3467 : PARALLEL_CATCH;
2134 :
2135 : // Set the cached NodalBCBase values in the Jacobian matrix
2136 3467 : _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{});
2137 3467 : }
2138 :
2139 : void
2140 717 : NonlinearSystemBase::getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs)
2141 : {
2142 717 : const Node & node = _mesh.nodeRef(node_id);
2143 717 : unsigned int s = number();
2144 717 : if (node.has_dofs(s))
2145 : {
2146 1472 : for (unsigned int v = 0; v < nVariables(); v++)
2147 1472 : for (unsigned int c = 0; c < node.n_comp(s, v); c++)
2148 717 : dofs.push_back(node.dof_number(s, v, c));
2149 : }
2150 717 : }
2151 :
2152 : void
2153 5567 : NonlinearSystemBase::findImplicitGeometricCouplingEntries(
2154 : GeometricSearchData & geom_search_data,
2155 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & graph)
2156 : {
2157 5567 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
2158 5567 : const auto & nearest_node_locators = geom_search_data._nearest_node_locators;
2159 7696 : for (const auto & it : nearest_node_locators)
2160 : {
2161 2129 : std::vector<dof_id_type> & secondary_nodes = it.second->_secondary_nodes;
2162 :
2163 10915 : for (const auto & secondary_node : secondary_nodes)
2164 : {
2165 8786 : std::set<dof_id_type> unique_secondary_indices;
2166 8786 : std::set<dof_id_type> unique_primary_indices;
2167 :
2168 8786 : auto node_to_elem_pair = node_to_elem_map.find(secondary_node);
2169 8786 : if (node_to_elem_pair != node_to_elem_map.end())
2170 : {
2171 8666 : const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
2172 :
2173 : // Get the dof indices from each elem connected to the node
2174 23314 : for (const auto & cur_elem : elems)
2175 : {
2176 14648 : std::vector<dof_id_type> dof_indices;
2177 14648 : dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2178 :
2179 73240 : for (const auto & dof : dof_indices)
2180 58592 : unique_secondary_indices.insert(dof);
2181 14648 : }
2182 : }
2183 :
2184 8786 : std::vector<dof_id_type> primary_nodes = it.second->_neighbor_nodes[secondary_node];
2185 :
2186 49576 : for (const auto & primary_node : primary_nodes)
2187 : {
2188 40790 : auto primary_node_to_elem_pair = node_to_elem_map.find(primary_node);
2189 : mooseAssert(primary_node_to_elem_pair != node_to_elem_map.end(),
2190 : "Missing entry in node to elem map");
2191 40790 : const std::vector<dof_id_type> & primary_node_elems = primary_node_to_elem_pair->second;
2192 :
2193 : // Get the dof indices from each elem connected to the node
2194 107078 : for (const auto & cur_elem : primary_node_elems)
2195 : {
2196 66288 : std::vector<dof_id_type> dof_indices;
2197 66288 : dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2198 :
2199 331440 : for (const auto & dof : dof_indices)
2200 265152 : unique_primary_indices.insert(dof);
2201 66288 : }
2202 : }
2203 :
2204 55414 : for (const auto & secondary_id : unique_secondary_indices)
2205 506684 : for (const auto & primary_id : unique_primary_indices)
2206 : {
2207 460056 : graph[secondary_id].push_back(primary_id);
2208 460056 : graph[primary_id].push_back(secondary_id);
2209 : }
2210 8786 : }
2211 : }
2212 :
2213 : // handle node-to-node constraints
2214 5567 : const auto & ncs = _constraints.getActiveNodalConstraints();
2215 5724 : for (const auto & nc : ncs)
2216 : {
2217 157 : std::vector<dof_id_type> primary_dofs;
2218 157 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
2219 314 : for (const auto & node_id : primary_node_ids)
2220 : {
2221 157 : Node * node = _mesh.queryNodePtr(node_id);
2222 157 : if (node && node->processor_id() == this->processor_id())
2223 : {
2224 142 : getNodeDofs(node_id, primary_dofs);
2225 : }
2226 : }
2227 :
2228 157 : _communicator.allgather(primary_dofs);
2229 :
2230 157 : std::vector<dof_id_type> secondary_dofs;
2231 157 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
2232 732 : for (const auto & node_id : secondary_node_ids)
2233 : {
2234 575 : Node * node = _mesh.queryNodePtr(node_id);
2235 575 : if (node && node->processor_id() == this->processor_id())
2236 : {
2237 575 : getNodeDofs(node_id, secondary_dofs);
2238 : }
2239 : }
2240 :
2241 157 : _communicator.allgather(secondary_dofs);
2242 :
2243 314 : for (const auto & primary_id : primary_dofs)
2244 870 : for (const auto & secondary_id : secondary_dofs)
2245 : {
2246 713 : graph[primary_id].push_back(secondary_id);
2247 713 : graph[secondary_id].push_back(primary_id);
2248 : }
2249 157 : }
2250 :
2251 : // Make every entry sorted and unique
2252 30319 : for (auto & it : graph)
2253 : {
2254 24752 : std::vector<dof_id_type> & row = it.second;
2255 24752 : std::sort(row.begin(), row.end());
2256 24752 : std::vector<dof_id_type>::iterator uit = std::unique(row.begin(), row.end());
2257 24752 : row.resize(uit - row.begin());
2258 : }
2259 5567 : }
2260 :
2261 : void
2262 4135 : NonlinearSystemBase::addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data)
2263 : {
2264 4135 : if (!hasMatrix(systemMatrixTag()))
2265 0 : mooseError("Need a system matrix ");
2266 :
2267 : // At this point, have no idea how to make
2268 : // this work with tag system
2269 4135 : auto & jacobian = getMatrix(systemMatrixTag());
2270 :
2271 4135 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
2272 :
2273 4135 : findImplicitGeometricCouplingEntries(geom_search_data, graph);
2274 :
2275 28010 : for (const auto & it : graph)
2276 : {
2277 23875 : dof_id_type dof = it.first;
2278 23875 : const auto & row = it.second;
2279 :
2280 281261 : for (const auto & coupled_dof : row)
2281 257386 : jacobian.add(dof, coupled_dof, 0);
2282 : }
2283 4135 : }
2284 :
2285 : void
2286 4763 : NonlinearSystemBase::constraintJacobians(const SparseMatrix<Number> & jacobian_to_view,
2287 : bool displaced)
2288 : {
2289 4763 : if (!hasMatrix(systemMatrixTag()))
2290 0 : mooseError("A system matrix is required");
2291 :
2292 4763 : auto & jacobian = getMatrix(systemMatrixTag());
2293 :
2294 4763 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2295 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2296 : MAT_NEW_NONZERO_ALLOCATION_ERR,
2297 : PETSC_FALSE));
2298 4763 : if (_fe_problem.ignoreZerosInJacobian())
2299 0 : LibmeshPetscCall(MatSetOption(
2300 : static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2301 :
2302 4763 : std::vector<numeric_index_type> zero_rows;
2303 :
2304 : if (displaced)
2305 : mooseAssert(_fe_problem.getDisplacedProblem(),
2306 : "If we're calling this method with displaced = true, then we better well have a "
2307 : "displaced problem");
2308 2192 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
2309 5859 : : static_cast<SubProblem &>(_fe_problem);
2310 4763 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
2311 :
2312 : bool constraints_applied;
2313 4763 : if (!_assemble_constraints_separately)
2314 4763 : constraints_applied = false;
2315 6829 : for (const auto & it : penetration_locators)
2316 : {
2317 2066 : if (_assemble_constraints_separately)
2318 : {
2319 : // Reset the constraint_applied flag before each new constraint, as they need to be
2320 : // assembled separately
2321 0 : constraints_applied = false;
2322 : }
2323 2066 : PenetrationLocator & pen_loc = *(it.second);
2324 :
2325 2066 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
2326 :
2327 2066 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
2328 2066 : BoundaryID primary_boundary = pen_loc._primary_boundary;
2329 :
2330 2066 : zero_rows.clear();
2331 2066 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
2332 : {
2333 : const auto & constraints =
2334 1128 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
2335 :
2336 5726 : for (const auto & secondary_node_num : secondary_nodes)
2337 : {
2338 4598 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
2339 :
2340 4598 : if (secondary_node.processor_id() == processor_id())
2341 : {
2342 4184 : if (pen_loc._penetration_info[secondary_node_num])
2343 : {
2344 4184 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
2345 :
2346 4184 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
2347 4184 : _fe_problem.reinitOffDiagScalars(0);
2348 :
2349 8368 : for (const auto & nfc : constraints)
2350 : {
2351 4184 : if (nfc->isExplicitConstraint())
2352 0 : continue;
2353 : // Return if this constraint does not correspond to the primary-secondary pair
2354 : // prepared by the outer loops.
2355 : // This continue statement is required when, e.g. one secondary surface constrains
2356 : // more than one primary surface.
2357 8368 : if (nfc->secondaryBoundary() != secondary_boundary ||
2358 4184 : nfc->primaryBoundary() != primary_boundary)
2359 0 : continue;
2360 :
2361 4184 : nfc->_jacobian = &jacobian_to_view;
2362 :
2363 4184 : if (nfc->shouldApply())
2364 : {
2365 4184 : constraints_applied = true;
2366 :
2367 4184 : nfc->prepareShapes(nfc->variable().number());
2368 4184 : nfc->prepareNeighborShapes(nfc->variable().number());
2369 :
2370 4184 : nfc->computeJacobian();
2371 :
2372 4184 : if (nfc->overwriteSecondaryJacobian())
2373 : {
2374 : // Add this variable's dof's row to be zeroed
2375 4184 : zero_rows.push_back(nfc->variable().nodalDofIndex());
2376 : }
2377 :
2378 4184 : std::vector<dof_id_type> secondary_dofs(1, nfc->variable().nodalDofIndex());
2379 :
2380 : // Assume that if the user is overwriting the secondary Jacobian, then they are
2381 : // supplying Jacobians that do not correspond to their other physics
2382 : // (e.g. Kernels), hence we should not apply a scalingFactor that is normally
2383 : // based on the order of their other physics (e.g. Kernels)
2384 : Real scaling_factor =
2385 4184 : nfc->overwriteSecondaryJacobian() ? 1. : nfc->variable().scalingFactor();
2386 :
2387 : // Cache the jacobian block for the secondary side
2388 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2389 4184 : nfc->_Kee,
2390 : secondary_dofs,
2391 4184 : nfc->_connected_dof_indices,
2392 : scaling_factor);
2393 :
2394 : // Cache Ken, Kne, Knn
2395 4184 : if (nfc->addCouplingEntriesToJacobian())
2396 : {
2397 : // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2398 : // factor when we're overwriting secondary stuff)
2399 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2400 4184 : nfc->_Ken,
2401 : secondary_dofs,
2402 4184 : nfc->primaryVariable().dofIndicesNeighbor(),
2403 : scaling_factor);
2404 :
2405 : // Use _connected_dof_indices to get all the correct columns
2406 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2407 4184 : nfc->_Kne,
2408 4184 : nfc->primaryVariable().dofIndicesNeighbor(),
2409 4184 : nfc->_connected_dof_indices,
2410 4184 : nfc->primaryVariable().scalingFactor());
2411 :
2412 : // We've handled Ken and Kne, finally handle Knn
2413 4184 : _fe_problem.cacheJacobianNeighbor(0);
2414 : }
2415 :
2416 : // Do the off-diagonals next
2417 4184 : const std::vector<MooseVariableFEBase *> coupled_vars = nfc->getCoupledMooseVars();
2418 8368 : for (const auto & jvar : coupled_vars)
2419 : {
2420 : // Only compute jacobians for nonlinear variables
2421 4184 : if (jvar->kind() != Moose::VAR_SOLVER)
2422 0 : continue;
2423 :
2424 : // Only compute Jacobian entries if this coupling is being used by the
2425 : // preconditioner
2426 4256 : if (nfc->variable().number() == jvar->number() ||
2427 144 : !_fe_problem.areCoupled(
2428 72 : nfc->variable().number(), jvar->number(), this->number()))
2429 4112 : continue;
2430 :
2431 : // Need to zero out the matrices first
2432 72 : _fe_problem.prepareAssembly(0);
2433 :
2434 72 : nfc->prepareShapes(nfc->variable().number());
2435 72 : nfc->prepareNeighborShapes(jvar->number());
2436 :
2437 72 : nfc->computeOffDiagJacobian(jvar->number());
2438 :
2439 : // Cache the jacobian block for the secondary side
2440 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2441 72 : nfc->_Kee,
2442 : secondary_dofs,
2443 72 : nfc->_connected_dof_indices,
2444 : scaling_factor);
2445 :
2446 : // Cache Ken, Kne, Knn
2447 72 : if (nfc->addCouplingEntriesToJacobian())
2448 : {
2449 : // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2450 : // factor when we're overwriting secondary stuff)
2451 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2452 72 : nfc->_Ken,
2453 : secondary_dofs,
2454 72 : jvar->dofIndicesNeighbor(),
2455 : scaling_factor);
2456 :
2457 : // Use _connected_dof_indices to get all the correct columns
2458 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2459 72 : nfc->_Kne,
2460 72 : nfc->variable().dofIndicesNeighbor(),
2461 72 : nfc->_connected_dof_indices,
2462 72 : nfc->variable().scalingFactor());
2463 :
2464 : // We've handled Ken and Kne, finally handle Knn
2465 72 : _fe_problem.cacheJacobianNeighbor(0);
2466 : }
2467 : }
2468 4184 : }
2469 : }
2470 : }
2471 : }
2472 : }
2473 : }
2474 2066 : if (_assemble_constraints_separately)
2475 : {
2476 : // See if constraints were applied anywhere
2477 0 : _communicator.max(constraints_applied);
2478 :
2479 0 : if (constraints_applied)
2480 : {
2481 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2482 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2483 : PETSC_TRUE));
2484 :
2485 0 : jacobian.close();
2486 0 : jacobian.zero_rows(zero_rows, 0.0);
2487 0 : jacobian.close();
2488 0 : _fe_problem.addCachedJacobian(0);
2489 0 : jacobian.close();
2490 : }
2491 : }
2492 : }
2493 4763 : if (!_assemble_constraints_separately)
2494 : {
2495 : // See if constraints were applied anywhere
2496 4763 : _communicator.max(constraints_applied);
2497 :
2498 4763 : if (constraints_applied)
2499 : {
2500 1054 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2501 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2502 : PETSC_TRUE));
2503 :
2504 1054 : jacobian.close();
2505 1054 : jacobian.zero_rows(zero_rows, 0.0);
2506 1054 : jacobian.close();
2507 1054 : _fe_problem.addCachedJacobian(0);
2508 1054 : jacobian.close();
2509 : }
2510 : }
2511 :
2512 4763 : THREAD_ID tid = 0;
2513 : // go over element-element constraint interface
2514 4763 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
2515 4763 : for (const auto & it : element_pair_locators)
2516 : {
2517 0 : ElementPairLocator & elem_pair_loc = *(it.second);
2518 :
2519 0 : if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
2520 : {
2521 : // ElemElemConstraint objects
2522 : const auto & _element_constraints =
2523 0 : _constraints.getActiveElemElemConstraints(it.first, displaced);
2524 :
2525 : // go over pair elements
2526 : const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
2527 0 : elem_pair_loc.getElemPairs();
2528 0 : for (const auto & pr : elem_pairs)
2529 : {
2530 0 : const Elem * elem1 = pr.first;
2531 0 : const Elem * elem2 = pr.second;
2532 :
2533 0 : if (elem1->processor_id() != processor_id())
2534 0 : continue;
2535 :
2536 0 : const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
2537 :
2538 : // for each element process constraints on the
2539 0 : for (const auto & ec : _element_constraints)
2540 : {
2541 0 : _fe_problem.setCurrentSubdomainID(elem1, tid);
2542 0 : subproblem.reinitElemPhys(elem1, info._elem1_constraint_q_point, tid);
2543 0 : _fe_problem.setNeighborSubdomainID(elem2, tid);
2544 0 : subproblem.reinitNeighborPhys(elem2, info._elem2_constraint_q_point, tid);
2545 :
2546 0 : ec->prepareShapes(ec->variable().number());
2547 0 : ec->prepareNeighborShapes(ec->variable().number());
2548 :
2549 0 : ec->reinit(info);
2550 0 : ec->computeJacobian();
2551 0 : _fe_problem.cacheJacobian(tid);
2552 0 : _fe_problem.cacheJacobianNeighbor(tid);
2553 : }
2554 0 : _fe_problem.addCachedJacobian(tid);
2555 : }
2556 : }
2557 : }
2558 :
2559 : // go over NodeElemConstraints
2560 4763 : std::set<dof_id_type> unique_secondary_node_ids;
2561 4763 : constraints_applied = false;
2562 20461 : for (const auto & secondary_id : _mesh.meshSubdomains())
2563 : {
2564 77282 : for (const auto & primary_id : _mesh.meshSubdomains())
2565 : {
2566 61584 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
2567 : {
2568 : const auto & constraints =
2569 180 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
2570 :
2571 : // get unique set of ids of all nodes on current block
2572 180 : unique_secondary_node_ids.clear();
2573 180 : const MeshBase & meshhelper = _mesh.getMesh();
2574 360 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
2575 8940 : meshhelper.active_subdomain_elements_end(secondary_id)))
2576 : {
2577 49140 : for (auto & n : elem->node_ref_range())
2578 40740 : unique_secondary_node_ids.insert(n.id());
2579 180 : }
2580 :
2581 14040 : for (auto secondary_node_id : unique_secondary_node_ids)
2582 : {
2583 13860 : const Node & secondary_node = _mesh.nodeRef(secondary_node_id);
2584 : // check if secondary node is on current processor
2585 13860 : if (secondary_node.processor_id() == processor_id())
2586 : {
2587 : // This reinits the variables that exist on the secondary node
2588 11088 : _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
2589 :
2590 : // This will set aside residual and jacobian space for the variables that have dofs
2591 : // on the secondary node
2592 11088 : _fe_problem.prepareAssembly(0);
2593 11088 : _fe_problem.reinitOffDiagScalars(0);
2594 :
2595 23584 : for (const auto & nec : constraints)
2596 : {
2597 12496 : if (nec->shouldApply())
2598 : {
2599 5680 : constraints_applied = true;
2600 :
2601 5680 : nec->_jacobian = &jacobian_to_view;
2602 5680 : nec->prepareShapes(nec->variable().number());
2603 5680 : nec->prepareNeighborShapes(nec->variable().number());
2604 :
2605 5680 : nec->computeJacobian();
2606 :
2607 5680 : if (nec->overwriteSecondaryJacobian())
2608 : {
2609 : // Add this variable's dof's row to be zeroed
2610 0 : zero_rows.push_back(nec->variable().nodalDofIndex());
2611 : }
2612 :
2613 5680 : std::vector<dof_id_type> secondary_dofs(1, nec->variable().nodalDofIndex());
2614 :
2615 : // Cache the jacobian block for the secondary side
2616 11360 : nec->addJacobian(_fe_problem.assembly(0, number()),
2617 5680 : nec->_Kee,
2618 : secondary_dofs,
2619 5680 : nec->_connected_dof_indices,
2620 5680 : nec->variable().scalingFactor());
2621 :
2622 : // Cache the jacobian block for the primary side
2623 11360 : nec->addJacobian(_fe_problem.assembly(0, number()),
2624 5680 : nec->_Kne,
2625 5680 : nec->primaryVariable().dofIndicesNeighbor(),
2626 5680 : nec->_connected_dof_indices,
2627 5680 : nec->primaryVariable().scalingFactor());
2628 :
2629 5680 : _fe_problem.cacheJacobian(0);
2630 5680 : _fe_problem.cacheJacobianNeighbor(0);
2631 :
2632 : // Do the off-diagonals next
2633 5680 : const std::vector<MooseVariableFEBase *> coupled_vars = nec->getCoupledMooseVars();
2634 11440 : for (const auto & jvar : coupled_vars)
2635 : {
2636 : // Only compute jacobians for nonlinear variables
2637 5760 : if (jvar->kind() != Moose::VAR_SOLVER)
2638 80 : continue;
2639 :
2640 : // Only compute Jacobian entries if this coupling is being used by the
2641 : // preconditioner
2642 6320 : if (nec->variable().number() == jvar->number() ||
2643 1280 : !_fe_problem.areCoupled(
2644 640 : nec->variable().number(), jvar->number(), this->number()))
2645 5040 : continue;
2646 :
2647 : // Need to zero out the matrices first
2648 640 : _fe_problem.prepareAssembly(0);
2649 :
2650 640 : nec->prepareShapes(nec->variable().number());
2651 640 : nec->prepareNeighborShapes(jvar->number());
2652 :
2653 640 : nec->computeOffDiagJacobian(jvar->number());
2654 :
2655 : // Cache the jacobian block for the secondary side
2656 1280 : nec->addJacobian(_fe_problem.assembly(0, number()),
2657 640 : nec->_Kee,
2658 : secondary_dofs,
2659 640 : nec->_connected_dof_indices,
2660 640 : nec->variable().scalingFactor());
2661 :
2662 : // Cache the jacobian block for the primary side
2663 1280 : nec->addJacobian(_fe_problem.assembly(0, number()),
2664 640 : nec->_Kne,
2665 640 : nec->variable().dofIndicesNeighbor(),
2666 640 : nec->_connected_dof_indices,
2667 640 : nec->variable().scalingFactor());
2668 :
2669 640 : _fe_problem.cacheJacobian(0);
2670 640 : _fe_problem.cacheJacobianNeighbor(0);
2671 : }
2672 5680 : }
2673 : }
2674 : }
2675 : }
2676 : }
2677 : }
2678 : }
2679 : // See if constraints were applied anywhere
2680 4763 : _communicator.max(constraints_applied);
2681 :
2682 4763 : if (constraints_applied)
2683 : {
2684 180 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2685 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2686 : PETSC_TRUE));
2687 :
2688 180 : jacobian.close();
2689 180 : jacobian.zero_rows(zero_rows, 0.0);
2690 180 : jacobian.close();
2691 180 : _fe_problem.addCachedJacobian(0);
2692 180 : jacobian.close();
2693 : }
2694 4763 : }
2695 :
2696 : void
2697 519689 : NonlinearSystemBase::computeScalarKernelsJacobians(const std::set<TagID> & tags)
2698 : {
2699 : MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
2700 :
2701 519689 : if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2702 510118 : scalar_kernel_warehouse = &_scalar_kernels;
2703 9571 : else if (tags.size() == 1)
2704 8319 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2705 : else
2706 1252 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagsObjectWarehouse(tags, 0));
2707 :
2708 : // Compute the diagonal block for scalar variables
2709 519689 : if (scalar_kernel_warehouse->hasActiveObjects())
2710 : {
2711 17554 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
2712 :
2713 17554 : _fe_problem.reinitScalars(/*tid=*/0);
2714 :
2715 17554 : _fe_problem.reinitOffDiagScalars(/*_tid*/ 0);
2716 :
2717 17554 : bool have_scalar_contributions = false;
2718 62904 : for (const auto & kernel : scalars)
2719 : {
2720 45350 : if (!kernel->computesJacobian())
2721 0 : continue;
2722 :
2723 45350 : kernel->reinit();
2724 45350 : const std::vector<dof_id_type> & dof_indices = kernel->variable().dofIndices();
2725 45350 : const DofMap & dof_map = kernel->variable().dofMap();
2726 45350 : const dof_id_type first_dof = dof_map.first_dof();
2727 45350 : const dof_id_type end_dof = dof_map.end_dof();
2728 51509 : for (dof_id_type dof : dof_indices)
2729 : {
2730 45392 : if (dof >= first_dof && dof < end_dof)
2731 : {
2732 39233 : kernel->computeJacobian();
2733 39233 : _fe_problem.addJacobianOffDiagScalar(kernel->variable().number());
2734 39233 : have_scalar_contributions = true;
2735 39233 : break;
2736 : }
2737 : }
2738 : }
2739 :
2740 17554 : if (have_scalar_contributions)
2741 15622 : _fe_problem.addJacobianScalar();
2742 : }
2743 519689 : }
2744 :
2745 : void
2746 535749 : NonlinearSystemBase::jacobianSetup()
2747 : {
2748 535749 : SolverSystem::jacobianSetup();
2749 :
2750 1120944 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2751 : {
2752 585195 : _kernels.jacobianSetup(tid);
2753 585195 : _nodal_kernels.jacobianSetup(tid);
2754 585195 : _dirac_kernels.jacobianSetup(tid);
2755 585195 : if (_doing_dg)
2756 2948 : _dg_kernels.jacobianSetup(tid);
2757 585195 : _interface_kernels.jacobianSetup(tid);
2758 585195 : _element_dampers.jacobianSetup(tid);
2759 585195 : _nodal_dampers.jacobianSetup(tid);
2760 585195 : _integrated_bcs.jacobianSetup(tid);
2761 : }
2762 535749 : _scalar_kernels.jacobianSetup();
2763 535749 : _constraints.jacobianSetup();
2764 535749 : _general_dampers.jacobianSetup();
2765 535749 : _nodal_bcs.jacobianSetup();
2766 :
2767 : // Avoid recursion
2768 535749 : if (this == &_fe_problem.currentNonlinearSystem())
2769 519689 : _fe_problem.jacobianSetup();
2770 535749 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
2771 535749 : }
2772 :
2773 : void
2774 519689 : NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
2775 : {
2776 1559067 : TIME_SECTION("computeJacobianInternal", 3);
2777 :
2778 519689 : _fe_problem.setCurrentNonlinearSystem(number());
2779 :
2780 : // Make matrix ready to use
2781 519689 : activateAllMatrixTags();
2782 :
2783 1551240 : for (auto tag : tags)
2784 : {
2785 1031551 : if (!hasMatrix(tag))
2786 510118 : continue;
2787 :
2788 521433 : auto & jacobian = getMatrix(tag);
2789 : // Necessary for speed
2790 521433 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2791 : {
2792 520147 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2793 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2794 : PETSC_TRUE));
2795 520147 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2796 13371 : LibmeshPetscCall(
2797 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2798 520147 : if (_fe_problem.ignoreZerosInJacobian())
2799 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2800 : MAT_IGNORE_ZERO_ENTRIES,
2801 : PETSC_TRUE));
2802 : }
2803 : }
2804 :
2805 519689 : jacobianSetup();
2806 :
2807 : // Jacobian contributions from UOs - for now this is used for ray tracing
2808 : // and ray kernels that contribute to the Jacobian (think line sources)
2809 519689 : std::vector<UserObject *> uos;
2810 519689 : _fe_problem.theWarehouse()
2811 1039378 : .query()
2812 519689 : .condition<AttribSystem>("UserObject")
2813 519689 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
2814 519689 : .queryInto(uos);
2815 519689 : for (auto & uo : uos)
2816 0 : uo->jacobianSetup();
2817 519689 : for (auto & uo : uos)
2818 : {
2819 0 : uo->initialize();
2820 0 : uo->execute();
2821 0 : uo->finalize();
2822 : }
2823 :
2824 : // reinit scalar variables
2825 1087451 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2826 567762 : _fe_problem.reinitScalars(tid);
2827 :
2828 : PARALLEL_TRY
2829 : {
2830 : // We would like to compute ScalarKernels, block NodalKernels, FVFluxKernels, and mortar objects
2831 : // up front because we want these included whether we are computing an ordinary Jacobian or a
2832 : // Jacobian for determining variable scaling factors
2833 519689 : computeScalarKernelsJacobians(tags);
2834 :
2835 : // Block restricted Nodal Kernels
2836 519689 : if (_nodal_kernels.hasActiveBlockObjects())
2837 : {
2838 3998 : ComputeNodalKernelJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
2839 3998 : const ConstNodeRange & range = _fe_problem.getCurrentAlgebraicNodeRange();
2840 3998 : Threads::parallel_reduce(range, cnkjt);
2841 :
2842 3998 : unsigned int n_threads = libMesh::n_threads();
2843 9262 : for (unsigned int i = 0; i < n_threads;
2844 : i++) // Add any cached jacobians that might be hanging around
2845 5264 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
2846 3998 : }
2847 :
2848 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
2849 519689 : if (_fe_problem.haveFV())
2850 : {
2851 : // the same loop works for both residual and jacobians because it keys
2852 : // off of FEProblem's _currently_computing_jacobian parameter
2853 : ComputeFVFluxJacobianThread<FVRange> fvj(
2854 37586 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
2855 37586 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
2856 37586 : Threads::parallel_reduce(faces, fvj);
2857 37586 : }
2858 519689 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
2859 519689 : displaced_problem && displaced_problem->haveFV())
2860 : {
2861 : ComputeFVFluxJacobianThread<FVRange> fvr(
2862 0 : _fe_problem, this->number(), tags, /*on_displaced=*/true);
2863 0 : FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
2864 0 : displaced_problem->mesh().ownedFaceInfoEnd());
2865 0 : Threads::parallel_reduce(faces, fvr);
2866 519689 : }
2867 :
2868 519689 : mortarConstraints(Moose::ComputeType::Jacobian, {}, tags);
2869 :
2870 : // Get our element range for looping over
2871 519689 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
2872 :
2873 519689 : if (_fe_problem.computingScalingJacobian())
2874 : {
2875 : // Only compute Jacobians corresponding to the diagonals of volumetric compute objects
2876 : // because this typically gives us a good representation of the physics. NodalBCs and
2877 : // Constraints can introduce dramatically different scales (often order unity).
2878 : // IntegratedBCs and/or InterfaceKernels may use penalty factors. DGKernels may be ok, but
2879 : // they are almost always used in conjunction with Kernels
2880 602 : ComputeJacobianForScalingThread cj(_fe_problem, tags);
2881 602 : Threads::parallel_reduce(elem_range, cj);
2882 602 : unsigned int n_threads = libMesh::n_threads();
2883 1263 : for (unsigned int i = 0; i < n_threads;
2884 : i++) // Add any Jacobian contributions still hanging around
2885 661 : _fe_problem.addCachedJacobian(i);
2886 :
2887 : // Check whether any exceptions were thrown and propagate this information for parallel
2888 : // consistency before
2889 : // 1) we do parallel communication when closing tagged matrices
2890 : // 2) early returning before reaching our PARALLEL_CATCH below
2891 602 : _fe_problem.checkExceptionAndStopSolve();
2892 :
2893 602 : closeTaggedMatrices(tags);
2894 :
2895 602 : return;
2896 602 : }
2897 :
2898 519087 : switch (_fe_problem.coupling())
2899 : {
2900 430269 : case Moose::COUPLING_DIAG:
2901 : {
2902 430269 : ComputeJacobianThread cj(_fe_problem, tags);
2903 430269 : Threads::parallel_reduce(elem_range, cj);
2904 :
2905 430261 : unsigned int n_threads = libMesh::n_threads();
2906 902236 : for (unsigned int i = 0; i < n_threads;
2907 : i++) // Add any Jacobian contributions still hanging around
2908 471975 : _fe_problem.addCachedJacobian(i);
2909 :
2910 : // Boundary restricted Nodal Kernels
2911 430261 : if (_nodal_kernels.hasActiveBoundaryObjects())
2912 : {
2913 46 : ComputeNodalKernelBCJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
2914 46 : const ConstBndNodeRange & bnd_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
2915 :
2916 46 : Threads::parallel_reduce(bnd_range, cnkjt);
2917 46 : unsigned int n_threads = libMesh::n_threads();
2918 96 : for (unsigned int i = 0; i < n_threads;
2919 : i++) // Add any cached jacobians that might be hanging around
2920 50 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
2921 46 : }
2922 430261 : }
2923 430261 : break;
2924 :
2925 88818 : default:
2926 : case Moose::COUPLING_CUSTOM:
2927 : {
2928 88818 : ComputeFullJacobianThread cj(_fe_problem, tags);
2929 88818 : Threads::parallel_reduce(elem_range, cj);
2930 88818 : unsigned int n_threads = libMesh::n_threads();
2931 :
2932 183931 : for (unsigned int i = 0; i < n_threads; i++)
2933 95117 : _fe_problem.addCachedJacobian(i);
2934 :
2935 : // Boundary restricted Nodal Kernels
2936 88814 : if (_nodal_kernels.hasActiveBoundaryObjects())
2937 : {
2938 10 : ComputeNodalKernelBCJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
2939 10 : const ConstBndNodeRange & bnd_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
2940 :
2941 10 : Threads::parallel_reduce(bnd_range, cnkjt);
2942 10 : unsigned int n_threads = libMesh::n_threads();
2943 21 : for (unsigned int i = 0; i < n_threads;
2944 : i++) // Add any cached jacobians that might be hanging around
2945 11 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
2946 10 : }
2947 88818 : }
2948 88814 : break;
2949 : }
2950 :
2951 519075 : computeDiracContributions(tags, true);
2952 :
2953 : static bool first = true;
2954 :
2955 : // This adds zeroes into geometric coupling entries to ensure they stay in the matrix
2956 1037526 : if ((_fe_problem.restoreOriginalNonzeroPattern() || first) &&
2957 518455 : _add_implicit_geometric_coupling_entries_to_jacobian)
2958 : {
2959 3051 : first = false;
2960 3051 : addImplicitGeometricCouplingEntries(_fe_problem.geomSearchData());
2961 :
2962 3051 : if (_fe_problem.getDisplacedProblem())
2963 1084 : addImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData());
2964 : }
2965 : }
2966 519071 : PARALLEL_CATCH;
2967 :
2968 : // Have no idea how to have constraints work
2969 : // with the tag system
2970 : PARALLEL_TRY
2971 : {
2972 : // Add in Jacobian contributions from other Constraints
2973 518969 : if (_fe_problem._has_constraints && tags.count(systemMatrixTag()))
2974 : {
2975 : // Some constraints need to be able to read values from the Jacobian, which requires that it
2976 : // be closed/assembled
2977 3667 : auto & system_matrix = getMatrix(systemMatrixTag());
2978 : #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
2979 : SparseMatrix<Number> * view_jac_ptr;
2980 3667 : std::unique_ptr<SparseMatrix<Number>> hash_copy;
2981 3667 : if (system_matrix.use_hash_table())
2982 : {
2983 2309 : hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
2984 2309 : view_jac_ptr = hash_copy.get();
2985 : }
2986 : else
2987 1358 : view_jac_ptr = &system_matrix;
2988 3667 : auto & jacobian_to_view = *view_jac_ptr;
2989 : #else
2990 : auto & jacobian_to_view = system_matrix;
2991 : #endif
2992 3667 : if (&jacobian_to_view == &system_matrix)
2993 1358 : system_matrix.close();
2994 :
2995 : // Nodal Constraints
2996 3667 : enforceNodalConstraintsJacobian();
2997 :
2998 : // Undisplaced Constraints
2999 3667 : constraintJacobians(jacobian_to_view, false);
3000 :
3001 : // Displaced Constraints
3002 3667 : if (_fe_problem.getDisplacedProblem())
3003 1096 : constraintJacobians(jacobian_to_view, true);
3004 3667 : }
3005 : }
3006 518969 : PARALLEL_CATCH;
3007 :
3008 : // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
3009 : // on boundary nodes
3010 518969 : if (_has_diag_save_in)
3011 185 : _fe_problem.getAuxiliarySystem().solution().close();
3012 :
3013 : PARALLEL_TRY
3014 : {
3015 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
3016 : // Select nodal kernels
3017 518969 : if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
3018 509408 : nbc_warehouse = &_nodal_bcs;
3019 9561 : else if (tags.size() == 1)
3020 8309 : nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
3021 : else
3022 1252 : nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0));
3023 :
3024 518969 : if (nbc_warehouse->hasActiveObjects())
3025 : {
3026 : // We may be switching from add to set. Moreover, we rely on a call to MatZeroRows to enforce
3027 : // the nodal boundary condition constraints which requires that the matrix be truly assembled
3028 : // as opposed to just flushed. Consequently we can't do the following despite any desire to
3029 : // keep our initial sparsity pattern honored (see https://gitlab.com/petsc/petsc/-/issues/852)
3030 : //
3031 : // flushTaggedMatrices(tags);
3032 427014 : closeTaggedMatrices(tags);
3033 :
3034 : // Cache the information about which BCs are coupled to which
3035 : // variables, so we don't have to figure it out for each node.
3036 427014 : std::map<std::string, std::set<unsigned int>> bc_involved_vars;
3037 427014 : const std::set<BoundaryID> & all_boundary_ids = _mesh.getBoundaryIDs();
3038 2147328 : for (const auto & bid : all_boundary_ids)
3039 : {
3040 : // Get reference to all the NodalBCs for this ID. This is only
3041 : // safe if there are NodalBCBases there to be gotten...
3042 1720314 : if (nbc_warehouse->hasActiveBoundaryObjects(bid))
3043 : {
3044 844916 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
3045 1774540 : for (const auto & bc : bcs)
3046 : {
3047 : const std::vector<MooseVariableFEBase *> & coupled_moose_vars =
3048 929624 : bc->getCoupledMooseVars();
3049 :
3050 : // Create the set of "involved" MOOSE nonlinear vars, which includes all coupled vars
3051 : // and the BC's own variable
3052 929624 : std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3053 987896 : for (const auto & coupled_var : coupled_moose_vars)
3054 58272 : if (coupled_var->kind() == Moose::VAR_SOLVER)
3055 56694 : var_set.insert(coupled_var->number());
3056 :
3057 929624 : var_set.insert(bc->variable().number());
3058 : }
3059 : }
3060 : }
3061 :
3062 : // reinit scalar variables again. This reinit does not re-fill any of the scalar variable
3063 : // solution arrays because that was done above. It only will reorder the derivative
3064 : // information for AD calculations to be suitable for NodalBC calculations
3065 892338 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3066 465324 : _fe_problem.reinitScalars(tid, true);
3067 :
3068 : // Get variable coupling list. We do all the NodalBCBase stuff on
3069 : // thread 0... The couplingEntries() data structure determines
3070 : // which variables are "coupled" as far as the preconditioner is
3071 : // concerned, not what variables a boundary condition specifically
3072 : // depends on.
3073 427014 : auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number());
3074 :
3075 : // Compute Jacobians for NodalBCBases
3076 427014 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3077 22372888 : for (const auto & bnode : bnd_nodes)
3078 : {
3079 21945874 : BoundaryID boundary_id = bnode->_bnd_id;
3080 21945874 : Node * node = bnode->_node;
3081 :
3082 32188157 : if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
3083 10242283 : node->processor_id() == processor_id())
3084 : {
3085 7987793 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3086 :
3087 7987793 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
3088 16991327 : for (const auto & bc : bcs)
3089 : {
3090 : // Get the set of involved MOOSE vars for this BC
3091 9003534 : std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3092 :
3093 : // Loop over all the variables whose Jacobian blocks are
3094 : // actually being computed, call computeOffDiagJacobian()
3095 : // for each one which is actually coupled (otherwise the
3096 : // value is zero.)
3097 22632240 : for (const auto & it : coupling_entries)
3098 : {
3099 13628706 : unsigned int ivar = it.first->number(), jvar = it.second->number();
3100 :
3101 : // We are only going to call computeOffDiagJacobian() if:
3102 : // 1.) the BC's variable is ivar
3103 : // 2.) jvar is "involved" with the BC (including jvar==ivar), and
3104 : // 3.) the BC should apply.
3105 13628706 : if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
3106 9092918 : bc->computeOffDiagJacobian(jvar);
3107 : }
3108 :
3109 9003534 : const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars();
3110 9003905 : for (const auto & jvariable : coupled_scalar_vars)
3111 371 : if (hasScalarVariable(jvariable->name()))
3112 371 : bc->computeOffDiagJacobianScalar(jvariable->number());
3113 : }
3114 : }
3115 : } // end loop over boundary nodes
3116 :
3117 : // Set the cached NodalBCBase values in the Jacobian matrix
3118 427014 : _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{});
3119 427014 : }
3120 : }
3121 518969 : PARALLEL_CATCH;
3122 :
3123 518969 : closeTaggedMatrices(tags);
3124 :
3125 : // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
3126 : // on boundary nodes
3127 518969 : if (_has_nodalbc_diag_save_in)
3128 23 : _fe_problem.getAuxiliarySystem().solution().close();
3129 :
3130 518969 : if (hasDiagSaveIn())
3131 185 : _fe_problem.getAuxiliarySystem().update();
3132 :
3133 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
3134 : // counters
3135 518969 : _app.solutionInvalidity().syncIteration();
3136 518969 : _app.solutionInvalidity().solutionInvalidAccumulation();
3137 520385 : }
3138 :
3139 : void
3140 0 : NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian)
3141 : {
3142 0 : _nl_matrix_tags.clear();
3143 :
3144 0 : auto & tags = _fe_problem.getMatrixTags();
3145 :
3146 0 : for (auto & tag : tags)
3147 0 : _nl_matrix_tags.insert(tag.second);
3148 :
3149 0 : computeJacobian(jacobian, _nl_matrix_tags);
3150 0 : }
3151 :
3152 : void
3153 0 : NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian, const std::set<TagID> & tags)
3154 : {
3155 0 : associateMatrixToTag(jacobian, systemMatrixTag());
3156 :
3157 0 : computeJacobianTags(tags);
3158 :
3159 0 : disassociateMatrixFromTag(jacobian, systemMatrixTag());
3160 0 : }
3161 :
3162 : void
3163 519689 : NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
3164 : {
3165 1559067 : TIME_SECTION("computeJacobianTags", 5);
3166 :
3167 519689 : FloatingPointExceptionGuard fpe_guard(_app);
3168 :
3169 : try
3170 : {
3171 519689 : computeJacobianInternal(tags);
3172 : }
3173 106 : catch (MooseException & e)
3174 : {
3175 : // The buck stops here, we have already handled the exception by
3176 : // calling stopSolve(), it is now up to PETSc to return a
3177 : // "diverged" reason during the next solve.
3178 102 : }
3179 519681 : }
3180 :
3181 : void
3182 228 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
3183 : {
3184 228 : _nl_matrix_tags.clear();
3185 :
3186 228 : auto & tags = _fe_problem.getMatrixTags();
3187 684 : for (auto & tag : tags)
3188 456 : _nl_matrix_tags.insert(tag.second);
3189 :
3190 228 : computeJacobianBlocks(blocks, _nl_matrix_tags);
3191 228 : }
3192 :
3193 : void
3194 588 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
3195 : const std::set<TagID> & tags)
3196 : {
3197 1764 : TIME_SECTION("computeJacobianBlocks", 3);
3198 588 : FloatingPointExceptionGuard fpe_guard(_app);
3199 :
3200 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3201 : {
3202 1512 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3203 :
3204 1512 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3205 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
3206 : PETSC_TRUE));
3207 1512 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
3208 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3209 : MAT_NEW_NONZERO_ALLOCATION_ERR,
3210 : PETSC_TRUE));
3211 :
3212 1512 : jacobian.zero();
3213 : }
3214 :
3215 1233 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3216 645 : _fe_problem.reinitScalars(tid);
3217 :
3218 : PARALLEL_TRY
3219 : {
3220 588 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
3221 588 : ComputeJacobianBlocksThread cjb(_fe_problem, blocks, tags);
3222 588 : Threads::parallel_reduce(elem_range, cjb);
3223 588 : }
3224 588 : PARALLEL_CATCH;
3225 :
3226 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3227 1512 : blocks[i]->_jacobian.close();
3228 :
3229 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3230 : {
3231 1512 : libMesh::System & precond_system = blocks[i]->_precond_system;
3232 1512 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3233 :
3234 1512 : unsigned int ivar = blocks[i]->_ivar;
3235 1512 : unsigned int jvar = blocks[i]->_jvar;
3236 :
3237 : // Dirichlet BCs
3238 1512 : std::vector<numeric_index_type> zero_rows;
3239 : PARALLEL_TRY
3240 : {
3241 1512 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3242 49712 : for (const auto & bnode : bnd_nodes)
3243 : {
3244 48200 : BoundaryID boundary_id = bnode->_bnd_id;
3245 48200 : Node * node = bnode->_node;
3246 :
3247 48200 : if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
3248 : {
3249 41544 : const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
3250 :
3251 41544 : if (node->processor_id() == processor_id())
3252 : {
3253 32468 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3254 :
3255 160190 : for (const auto & bc : bcs)
3256 127722 : if (bc->variable().number() == ivar && bc->shouldApply())
3257 : {
3258 : // The first zero is for the variable number... there is only one variable in
3259 : // each mini-system The second zero only works with Lagrange elements!
3260 47568 : zero_rows.push_back(node->dof_number(precond_system.number(), 0, 0));
3261 : }
3262 : }
3263 : }
3264 : }
3265 : }
3266 1512 : PARALLEL_CATCH;
3267 :
3268 1512 : jacobian.close();
3269 :
3270 : // This zeroes the rows corresponding to Dirichlet BCs and puts a 1.0 on the diagonal
3271 1512 : if (ivar == jvar)
3272 1396 : jacobian.zero_rows(zero_rows, 1.0);
3273 : else
3274 116 : jacobian.zero_rows(zero_rows, 0.0);
3275 :
3276 1512 : jacobian.close();
3277 1512 : }
3278 588 : }
3279 :
3280 : void
3281 365967 : NonlinearSystemBase::updateActive(THREAD_ID tid)
3282 : {
3283 365967 : _element_dampers.updateActive(tid);
3284 365967 : _nodal_dampers.updateActive(tid);
3285 365967 : _integrated_bcs.updateActive(tid);
3286 365967 : _dg_kernels.updateActive(tid);
3287 365967 : _interface_kernels.updateActive(tid);
3288 365967 : _dirac_kernels.updateActive(tid);
3289 365967 : _kernels.updateActive(tid);
3290 365967 : _nodal_kernels.updateActive(tid);
3291 365967 : if (tid == 0)
3292 : {
3293 336220 : _general_dampers.updateActive();
3294 336220 : _nodal_bcs.updateActive();
3295 336220 : _preset_nodal_bcs.updateActive();
3296 336220 : _ad_preset_nodal_bcs.updateActive();
3297 336220 : _constraints.updateActive();
3298 336220 : _scalar_kernels.updateActive();
3299 : }
3300 365967 : }
3301 :
3302 : Real
3303 1420 : NonlinearSystemBase::computeDamping(const NumericVector<Number> & solution,
3304 : const NumericVector<Number> & update)
3305 : {
3306 : // Default to no damping
3307 1420 : Real damping = 1.0;
3308 1420 : bool has_active_dampers = false;
3309 :
3310 : try
3311 : {
3312 1420 : if (_element_dampers.hasActiveObjects())
3313 : {
3314 : PARALLEL_TRY
3315 : {
3316 2450 : TIME_SECTION("computeDampers", 3, "Computing Dampers");
3317 490 : has_active_dampers = true;
3318 490 : *_increment_vec = update;
3319 490 : ComputeElemDampingThread cid(_fe_problem, *this);
3320 490 : Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicElementRange(), cid);
3321 490 : damping = std::min(cid.damping(), damping);
3322 490 : }
3323 490 : PARALLEL_CATCH;
3324 : }
3325 :
3326 1370 : if (_nodal_dampers.hasActiveObjects())
3327 : {
3328 : PARALLEL_TRY
3329 : {
3330 2620 : TIME_SECTION("computeDamping::element", 3, "Computing Element Damping");
3331 :
3332 524 : has_active_dampers = true;
3333 524 : *_increment_vec = update;
3334 524 : ComputeNodalDampingThread cndt(_fe_problem, *this);
3335 524 : Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicNodeRange(), cndt);
3336 524 : damping = std::min(cndt.damping(), damping);
3337 524 : }
3338 524 : PARALLEL_CATCH;
3339 : }
3340 :
3341 1314 : if (_general_dampers.hasActiveObjects())
3342 : {
3343 : PARALLEL_TRY
3344 : {
3345 2020 : TIME_SECTION("computeDamping::general", 3, "Computing General Damping");
3346 :
3347 404 : has_active_dampers = true;
3348 404 : const auto & gdampers = _general_dampers.getActiveObjects();
3349 808 : for (const auto & damper : gdampers)
3350 : {
3351 404 : Real gd_damping = damper->computeDamping(solution, update);
3352 : try
3353 : {
3354 404 : damper->checkMinDamping(gd_damping);
3355 : }
3356 8 : catch (MooseException & e)
3357 : {
3358 16 : _fe_problem.setException(e.what());
3359 8 : }
3360 404 : damping = std::min(gd_damping, damping);
3361 : }
3362 404 : }
3363 404 : PARALLEL_CATCH;
3364 : }
3365 : }
3366 114 : catch (MooseException & e)
3367 : {
3368 : // The buck stops here, we have already handled the exception by
3369 : // calling stopSolve(), it is now up to PETSc to return a
3370 : // "diverged" reason during the next solve.
3371 114 : }
3372 :
3373 1420 : _communicator.min(damping);
3374 :
3375 1420 : if (has_active_dampers && damping < 1.0)
3376 995 : _console << " Damping factor: " << damping << std::endl;
3377 :
3378 1420 : return damping;
3379 : }
3380 :
3381 : void
3382 3812352 : NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
3383 : {
3384 3812352 : _fe_problem.clearDiracInfo();
3385 :
3386 3812352 : std::set<const Elem *> dirac_elements;
3387 :
3388 3812352 : if (_dirac_kernels.hasActiveObjects())
3389 : {
3390 192085 : TIME_SECTION("computeDirac", 3, "Computing DiracKernels");
3391 :
3392 : // TODO: Need a threading fix... but it's complicated!
3393 79882 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3394 : {
3395 41481 : const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
3396 96376 : for (const auto & dkernel : dkernels)
3397 : {
3398 54911 : dkernel->clearPoints();
3399 54911 : dkernel->addPoints();
3400 : }
3401 : }
3402 :
3403 38401 : ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
3404 :
3405 38401 : _fe_problem.getDiracElements(dirac_elements);
3406 :
3407 38401 : DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
3408 : // TODO: Make Dirac work thread!
3409 : // Threads::parallel_reduce(range, cd);
3410 :
3411 38401 : cd(range);
3412 38401 : }
3413 3812336 : }
3414 :
3415 : NumericVector<Number> &
3416 0 : NonlinearSystemBase::residualCopy()
3417 : {
3418 0 : if (!_residual_copy.get())
3419 0 : _residual_copy = NumericVector<Number>::build(_communicator);
3420 :
3421 0 : return *_residual_copy;
3422 : }
3423 :
3424 : NumericVector<Number> &
3425 360 : NonlinearSystemBase::residualGhosted()
3426 : {
3427 360 : _need_residual_ghosted = true;
3428 360 : if (!_residual_ghosted)
3429 : {
3430 : // The first time we realize we need a ghosted residual vector,
3431 : // we add it.
3432 576 : _residual_ghosted = &addVector("residual_ghosted", false, GHOSTED);
3433 :
3434 : // If we've already realized we need time and/or non-time
3435 : // residual vectors, but we haven't yet realized they need to be
3436 : // ghosted, fix that now.
3437 : //
3438 : // If an application changes its mind, the libMesh API lets us
3439 : // change the vector.
3440 288 : if (_Re_time)
3441 : {
3442 50 : const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
3443 50 : _Re_time = &system().add_vector(vector_name, false, GHOSTED);
3444 50 : }
3445 288 : if (_Re_non_time)
3446 : {
3447 288 : const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
3448 288 : _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
3449 288 : }
3450 : }
3451 360 : return *_residual_ghosted;
3452 : }
3453 :
3454 : void
3455 68579 : NonlinearSystemBase::augmentSparsity(SparsityPattern::Graph & sparsity,
3456 : std::vector<dof_id_type> & n_nz,
3457 : std::vector<dof_id_type> & n_oz)
3458 : {
3459 68579 : if (_add_implicit_geometric_coupling_entries_to_jacobian)
3460 : {
3461 1231 : _fe_problem.updateGeomSearch();
3462 :
3463 1231 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
3464 :
3465 1231 : findImplicitGeometricCouplingEntries(_fe_problem.geomSearchData(), graph);
3466 :
3467 1231 : if (_fe_problem.getDisplacedProblem())
3468 201 : findImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData(),
3469 : graph);
3470 :
3471 1231 : const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
3472 1231 : const dof_id_type end_dof_on_proc = dofMap().end_dof(processor_id());
3473 :
3474 : // The total number of dofs on and off processor
3475 1231 : const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
3476 1231 : const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
3477 :
3478 2088 : for (const auto & git : graph)
3479 : {
3480 857 : dof_id_type dof = git.first;
3481 857 : dof_id_type local_dof = dof - first_dof_on_proc;
3482 :
3483 857 : if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
3484 126 : continue;
3485 :
3486 731 : const auto & row = git.second;
3487 :
3488 731 : SparsityPattern::Row & sparsity_row = sparsity[local_dof];
3489 :
3490 731 : unsigned int original_row_length = sparsity_row.size();
3491 :
3492 731 : sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
3493 :
3494 1462 : SparsityPattern::sort_row(
3495 731 : sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
3496 :
3497 : // Fix up nonzero arrays
3498 4807 : for (const auto & coupled_dof : row)
3499 : {
3500 4076 : if (coupled_dof < first_dof_on_proc || coupled_dof >= end_dof_on_proc)
3501 : {
3502 1224 : if (n_oz[local_dof] < n_dofs_not_on_proc)
3503 612 : n_oz[local_dof]++;
3504 : }
3505 : else
3506 : {
3507 3464 : if (n_nz[local_dof] < n_dofs_on_proc)
3508 3464 : n_nz[local_dof]++;
3509 : }
3510 : }
3511 : }
3512 1231 : }
3513 68579 : }
3514 :
3515 : void
3516 0 : NonlinearSystemBase::setSolutionUDot(const NumericVector<Number> & u_dot)
3517 : {
3518 0 : *_u_dot = u_dot;
3519 0 : }
3520 :
3521 : void
3522 0 : NonlinearSystemBase::setSolutionUDotDot(const NumericVector<Number> & u_dotdot)
3523 : {
3524 0 : *_u_dotdot = u_dotdot;
3525 0 : }
3526 :
3527 : void
3528 0 : NonlinearSystemBase::setSolutionUDotOld(const NumericVector<Number> & u_dot_old)
3529 : {
3530 0 : *_u_dot_old = u_dot_old;
3531 0 : }
3532 :
3533 : void
3534 0 : NonlinearSystemBase::setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old)
3535 : {
3536 0 : *_u_dotdot_old = u_dotdot_old;
3537 0 : }
3538 :
3539 : void
3540 14633 : NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
3541 : {
3542 14633 : if (_preconditioner.get() != nullptr)
3543 7 : mooseError("More than one active Preconditioner detected");
3544 :
3545 14626 : _preconditioner = pc;
3546 14626 : }
3547 :
3548 : MoosePreconditioner const *
3549 54384 : NonlinearSystemBase::getPreconditioner() const
3550 : {
3551 54384 : return _preconditioner.get();
3552 : }
3553 :
3554 : void
3555 322 : NonlinearSystemBase::setupDampers()
3556 : {
3557 322 : _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
3558 322 : }
3559 :
3560 : void
3561 1753 : NonlinearSystemBase::reinitIncrementAtQpsForDampers(THREAD_ID /*tid*/,
3562 : const std::set<MooseVariable *> & damped_vars)
3563 : {
3564 3529 : for (const auto & var : damped_vars)
3565 1776 : var->computeIncrementAtQps(*_increment_vec);
3566 1753 : }
3567 :
3568 : void
3569 7441 : NonlinearSystemBase::reinitIncrementAtNodeForDampers(THREAD_ID /*tid*/,
3570 : const std::set<MooseVariable *> & damped_vars)
3571 : {
3572 14882 : for (const auto & var : damped_vars)
3573 7441 : var->computeIncrementAtNode(*_increment_vec);
3574 7441 : }
3575 :
3576 : void
3577 42945 : NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
3578 : {
3579 : // Obtain all blocks and variables covered by all kernels
3580 42945 : std::set<SubdomainID> input_subdomains;
3581 42945 : std::set<std::string> kernel_variables;
3582 :
3583 42945 : bool global_kernels_exist = false;
3584 42945 : global_kernels_exist |= _scalar_kernels.hasActiveObjects();
3585 42945 : global_kernels_exist |= _nodal_kernels.hasActiveObjects();
3586 :
3587 42945 : _kernels.subdomainsCovered(input_subdomains, kernel_variables);
3588 42945 : _dg_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3589 42945 : _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3590 42945 : _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3591 42945 : _constraints.subdomainsCovered(input_subdomains, kernel_variables);
3592 :
3593 42945 : if (_fe_problem.haveFV())
3594 : {
3595 3525 : std::vector<FVElementalKernel *> fv_elemental_kernels;
3596 3525 : _fe_problem.theWarehouse()
3597 7050 : .query()
3598 3525 : .template condition<AttribSystem>("FVElementalKernel")
3599 3525 : .queryInto(fv_elemental_kernels);
3600 :
3601 7533 : for (auto fv_kernel : fv_elemental_kernels)
3602 : {
3603 4008 : if (fv_kernel->blockRestricted())
3604 1972 : for (auto block_id : fv_kernel->blockIDs())
3605 1103 : input_subdomains.insert(block_id);
3606 : else
3607 3139 : global_kernels_exist = true;
3608 4008 : kernel_variables.insert(fv_kernel->variable().name());
3609 :
3610 : // Check for lagrange multiplier
3611 4008 : if (dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel))
3612 268 : kernel_variables.insert(dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel)
3613 268 : ->lambdaVariable()
3614 134 : .name());
3615 : }
3616 :
3617 3525 : std::vector<FVFluxKernel *> fv_flux_kernels;
3618 3525 : _fe_problem.theWarehouse()
3619 7050 : .query()
3620 3525 : .template condition<AttribSystem>("FVFluxKernel")
3621 3525 : .queryInto(fv_flux_kernels);
3622 :
3623 8792 : for (auto fv_kernel : fv_flux_kernels)
3624 : {
3625 5267 : if (fv_kernel->blockRestricted())
3626 2350 : for (auto block_id : fv_kernel->blockIDs())
3627 1250 : input_subdomains.insert(block_id);
3628 : else
3629 4167 : global_kernels_exist = true;
3630 5267 : kernel_variables.insert(fv_kernel->variable().name());
3631 : }
3632 :
3633 3525 : std::vector<FVInterfaceKernel *> fv_interface_kernels;
3634 3525 : _fe_problem.theWarehouse()
3635 7050 : .query()
3636 3525 : .template condition<AttribSystem>("FVInterfaceKernel")
3637 3525 : .queryInto(fv_interface_kernels);
3638 :
3639 3830 : for (auto fvik : fv_interface_kernels)
3640 305 : if (auto scalar_fvik = dynamic_cast<FVScalarLagrangeMultiplierInterface *>(fvik))
3641 14 : kernel_variables.insert(scalar_fvik->lambdaVariable().name());
3642 :
3643 3525 : std::vector<FVFluxBC *> fv_flux_bcs;
3644 3525 : _fe_problem.theWarehouse()
3645 7050 : .query()
3646 3525 : .template condition<AttribSystem>("FVFluxBC")
3647 3525 : .queryInto(fv_flux_bcs);
3648 :
3649 5555 : for (auto fvbc : fv_flux_bcs)
3650 2030 : if (auto scalar_fvbc = dynamic_cast<FVBoundaryScalarLagrangeMultiplierConstraint *>(fvbc))
3651 14 : kernel_variables.insert(scalar_fvbc->lambdaVariable().name());
3652 3525 : }
3653 :
3654 : // Check kernel coverage of subdomains (blocks) in your mesh
3655 42945 : if (!global_kernels_exist)
3656 : {
3657 38942 : std::set<SubdomainID> difference;
3658 38942 : std::set_difference(mesh_subdomains.begin(),
3659 : mesh_subdomains.end(),
3660 : input_subdomains.begin(),
3661 : input_subdomains.end(),
3662 : std::inserter(difference, difference.end()));
3663 :
3664 : // there supposed to be no kernels on this lower-dimensional subdomain
3665 39180 : for (const auto & id : _mesh.interiorLowerDBlocks())
3666 238 : difference.erase(id);
3667 39180 : for (const auto & id : _mesh.boundaryLowerDBlocks())
3668 238 : difference.erase(id);
3669 :
3670 38942 : if (!difference.empty())
3671 : {
3672 : std::vector<SubdomainID> difference_vec =
3673 8 : std::vector<SubdomainID>(difference.begin(), difference.end());
3674 8 : std::vector<SubdomainName> difference_names = _mesh.getSubdomainNames(difference_vec);
3675 8 : std::stringstream missing_block_names;
3676 8 : std::copy(difference_names.begin(),
3677 : difference_names.end(),
3678 8 : std::ostream_iterator<std::string>(missing_block_names, " "));
3679 8 : std::stringstream missing_block_ids;
3680 8 : std::copy(difference.begin(),
3681 : difference.end(),
3682 8 : std::ostream_iterator<unsigned int>(missing_block_ids, " "));
3683 :
3684 8 : mooseError("Each subdomain must contain at least one Kernel.\nThe following block(s) lack an "
3685 8 : "active kernel: " +
3686 8 : missing_block_names.str(),
3687 : " (ids: ",
3688 8 : missing_block_ids.str(),
3689 : ")");
3690 0 : }
3691 38934 : }
3692 :
3693 : // Check kernel use of variables
3694 42937 : std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
3695 :
3696 42937 : std::set<VariableName> difference;
3697 42937 : std::set_difference(variables.begin(),
3698 : variables.end(),
3699 : kernel_variables.begin(),
3700 : kernel_variables.end(),
3701 : std::inserter(difference, difference.end()));
3702 :
3703 : // skip checks for varaibles defined on lower-dimensional subdomain
3704 42937 : std::set<VariableName> vars(difference);
3705 43300 : for (auto & var_name : vars)
3706 : {
3707 363 : auto blks = getSubdomainsForVar(var_name);
3708 744 : for (const auto & id : blks)
3709 381 : if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
3710 381 : difference.erase(var_name);
3711 363 : }
3712 :
3713 42937 : if (!difference.empty())
3714 : {
3715 8 : std::stringstream missing_kernel_vars;
3716 8 : std::copy(difference.begin(),
3717 : difference.end(),
3718 8 : std::ostream_iterator<std::string>(missing_kernel_vars, " "));
3719 8 : mooseError("Each variable must be referenced by at least one active Kernel.\nThe following "
3720 8 : "variable(s) lack an active kernel: " +
3721 8 : missing_kernel_vars.str());
3722 0 : }
3723 42929 : }
3724 :
3725 : bool
3726 28854 : NonlinearSystemBase::containsTimeKernel()
3727 : {
3728 28854 : auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3729 :
3730 28854 : return time_kernels.hasActiveObjects();
3731 : }
3732 :
3733 : std::vector<std::string>
3734 0 : NonlinearSystemBase::timeKernelVariableNames()
3735 : {
3736 0 : std::vector<std::string> variable_names;
3737 0 : const auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3738 0 : if (time_kernels.hasActiveObjects())
3739 0 : for (const auto & kernel : time_kernels.getObjects())
3740 0 : variable_names.push_back(kernel->variable().name());
3741 :
3742 0 : return variable_names;
3743 0 : }
3744 :
3745 : bool
3746 28897 : NonlinearSystemBase::needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3747 : {
3748 28897 : return _integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid);
3749 : }
3750 :
3751 : bool
3752 1478 : NonlinearSystemBase::needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3753 : {
3754 1478 : return _interface_kernels.hasActiveBoundaryObjects(bnd_id, tid);
3755 : }
3756 :
3757 : bool
3758 12872 : NonlinearSystemBase::needSubdomainMaterialOnSide(SubdomainID /*subdomain_id*/,
3759 : THREAD_ID /*tid*/) const
3760 : {
3761 12872 : return _doing_dg;
3762 : }
3763 :
3764 : bool
3765 0 : NonlinearSystemBase::doingDG() const
3766 : {
3767 0 : return _doing_dg;
3768 : }
3769 :
3770 : void
3771 556 : NonlinearSystemBase::setPreviousNewtonSolution(const NumericVector<Number> & soln)
3772 : {
3773 556 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
3774 556 : getVector(Moose::PREVIOUS_NL_SOLUTION_TAG) = soln;
3775 556 : }
3776 :
3777 : void
3778 3816433 : NonlinearSystemBase::mortarConstraints(const Moose::ComputeType compute_type,
3779 : const std::set<TagID> & vector_tags,
3780 : const std::set<TagID> & matrix_tags)
3781 : {
3782 : parallel_object_only();
3783 :
3784 : try
3785 : {
3786 3825552 : for (auto & map_pr : _undisplaced_mortar_functors)
3787 9119 : map_pr.second(compute_type, vector_tags, matrix_tags);
3788 :
3789 3821121 : for (auto & map_pr : _displaced_mortar_functors)
3790 4688 : map_pr.second(compute_type, vector_tags, matrix_tags);
3791 : }
3792 0 : catch (MetaPhysicL::LogicError &)
3793 : {
3794 0 : mooseError(
3795 : "We caught a MetaPhysicL error in NonlinearSystemBase::mortarConstraints. This is very "
3796 : "likely due to AD not having a sufficiently large derivative container size. Please run "
3797 : "MOOSE configure with the '--with-derivative-size=<n>' option");
3798 0 : }
3799 3816433 : }
3800 :
3801 : void
3802 439 : NonlinearSystemBase::setupScalingData()
3803 : {
3804 439 : if (_auto_scaling_initd)
3805 0 : return;
3806 :
3807 : // Want the libMesh count of variables, not MOOSE, e.g. I don't care about array variable counts
3808 439 : const auto n_vars = system().n_vars();
3809 :
3810 439 : if (_scaling_group_variables.empty())
3811 : {
3812 427 : _var_to_group_var.reserve(n_vars);
3813 427 : _num_scaling_groups = n_vars;
3814 :
3815 1128 : for (const auto var_number : make_range(n_vars))
3816 701 : _var_to_group_var.emplace(var_number, var_number);
3817 : }
3818 : else
3819 : {
3820 12 : std::set<unsigned int> var_numbers, var_numbers_covered, var_numbers_not_covered;
3821 48 : for (const auto var_number : make_range(n_vars))
3822 36 : var_numbers.insert(var_number);
3823 :
3824 12 : _num_scaling_groups = _scaling_group_variables.size();
3825 :
3826 24 : for (const auto group_index : index_range(_scaling_group_variables))
3827 48 : for (const auto & var_name : _scaling_group_variables[group_index])
3828 : {
3829 36 : if (!hasVariable(var_name) && !hasScalarVariable(var_name))
3830 0 : mooseError("'",
3831 : var_name,
3832 : "', provided to the 'scaling_group_variables' parameter, does not exist in "
3833 : "the nonlinear system.");
3834 :
3835 : const MooseVariableBase & var =
3836 36 : hasVariable(var_name)
3837 24 : ? static_cast<MooseVariableBase &>(getVariable(0, var_name))
3838 60 : : static_cast<MooseVariableBase &>(getScalarVariable(0, var_name));
3839 36 : auto map_pair = _var_to_group_var.emplace(var.number(), group_index);
3840 36 : if (!map_pair.second)
3841 0 : mooseError("Variable ", var_name, " is contained in multiple scaling grouplings");
3842 36 : var_numbers_covered.insert(var.number());
3843 : }
3844 :
3845 12 : std::set_difference(var_numbers.begin(),
3846 : var_numbers.end(),
3847 : var_numbers_covered.begin(),
3848 : var_numbers_covered.end(),
3849 : std::inserter(var_numbers_not_covered, var_numbers_not_covered.begin()));
3850 :
3851 12 : _num_scaling_groups = _scaling_group_variables.size() + var_numbers_not_covered.size();
3852 :
3853 12 : auto index = static_cast<unsigned int>(_scaling_group_variables.size());
3854 12 : for (auto var_number : var_numbers_not_covered)
3855 0 : _var_to_group_var.emplace(var_number, index++);
3856 12 : }
3857 :
3858 439 : _variable_autoscaled.resize(n_vars, true);
3859 439 : const auto & number_to_var_map = _vars[0].numberToVariableMap();
3860 :
3861 439 : if (_ignore_variables_for_autoscaling.size())
3862 50 : for (const auto i : index_range(_variable_autoscaled))
3863 40 : if (std::find(_ignore_variables_for_autoscaling.begin(),
3864 : _ignore_variables_for_autoscaling.end(),
3865 80 : libmesh_map_find(number_to_var_map, i)->name()) !=
3866 80 : _ignore_variables_for_autoscaling.end())
3867 20 : _variable_autoscaled[i] = false;
3868 :
3869 439 : _auto_scaling_initd = true;
3870 : }
3871 :
3872 : bool
3873 1112 : NonlinearSystemBase::computeScaling()
3874 : {
3875 1112 : if (_compute_scaling_once && _computed_scaling)
3876 460 : return true;
3877 :
3878 652 : _console << "\nPerforming automatic scaling calculation\n" << std::endl;
3879 :
3880 3260 : TIME_SECTION("computeScaling", 3, "Computing Automatic Scaling");
3881 :
3882 : // It's funny but we need to assemble our vector of scaling factors here otherwise we will be
3883 : // applying scaling factors of 0 during Assembly of our scaling Jacobian
3884 652 : assembleScalingVector();
3885 :
3886 : // container for repeated access of element global dof indices
3887 652 : std::vector<dof_id_type> dof_indices;
3888 :
3889 652 : if (!_auto_scaling_initd)
3890 439 : setupScalingData();
3891 :
3892 652 : std::vector<Real> inverse_scaling_factors(_num_scaling_groups, 0);
3893 652 : std::vector<Real> resid_inverse_scaling_factors(_num_scaling_groups, 0);
3894 652 : std::vector<Real> jac_inverse_scaling_factors(_num_scaling_groups, 0);
3895 652 : auto & dof_map = dofMap();
3896 :
3897 : // what types of scaling do we want?
3898 652 : bool jac_scaling = _resid_vs_jac_scaling_param < 1. - TOLERANCE;
3899 652 : bool resid_scaling = _resid_vs_jac_scaling_param > TOLERANCE;
3900 :
3901 652 : const NumericVector<Number> & scaling_residual = RHS();
3902 :
3903 652 : if (jac_scaling)
3904 : {
3905 : // if (!_auto_scaling_initd)
3906 : // We need to reinit this when the number of dofs changes
3907 : // but there is no good way to track that
3908 : // In theory, it is the job of libmesh system to track this,
3909 : // but this special matrix is not owned by libMesh system
3910 : // Let us reinit eveytime since it is not expensive
3911 : {
3912 602 : auto init_vector = NumericVector<Number>::build(this->comm());
3913 602 : init_vector->init(system().n_dofs(), system().n_local_dofs(), /*fast=*/false, PARALLEL);
3914 :
3915 602 : _scaling_matrix->clear();
3916 602 : _scaling_matrix->init(*init_vector);
3917 602 : }
3918 :
3919 602 : _fe_problem.computingScalingJacobian(true);
3920 : // Dispatch to derived classes to ensure that we use the correct matrix tag
3921 602 : computeScalingJacobian();
3922 602 : _fe_problem.computingScalingJacobian(false);
3923 : }
3924 :
3925 652 : if (resid_scaling)
3926 : {
3927 50 : _fe_problem.computingScalingResidual(true);
3928 50 : _fe_problem.computingNonlinearResid(true);
3929 : // Dispatch to derived classes to ensure that we use the correct vector tag
3930 50 : computeScalingResidual();
3931 50 : _fe_problem.computingNonlinearResid(false);
3932 50 : _fe_problem.computingScalingResidual(false);
3933 : }
3934 :
3935 : // Did something bad happen during residual/Jacobian scaling computation?
3936 652 : if (_fe_problem.getFailNextNonlinearConvergenceCheck())
3937 0 : return false;
3938 :
3939 228330 : auto examine_dof_indices = [this,
3940 : jac_scaling,
3941 : resid_scaling,
3942 : &dof_map,
3943 : &jac_inverse_scaling_factors,
3944 : &resid_inverse_scaling_factors,
3945 : &scaling_residual](const auto & dof_indices, const auto var_number)
3946 : {
3947 1136309 : for (auto dof_index : dof_indices)
3948 907979 : if (dof_map.local_index(dof_index))
3949 : {
3950 903220 : if (jac_scaling)
3951 : {
3952 : // For now we will use the diagonal for determining scaling
3953 899103 : auto mat_value = (*_scaling_matrix)(dof_index, dof_index);
3954 899103 : auto & factor = jac_inverse_scaling_factors[_var_to_group_var[var_number]];
3955 899103 : factor = std::max(factor, std::abs(mat_value));
3956 : }
3957 903220 : if (resid_scaling)
3958 : {
3959 4117 : auto vec_value = scaling_residual(dof_index);
3960 4117 : auto & factor = resid_inverse_scaling_factors[_var_to_group_var[var_number]];
3961 4117 : factor = std::max(factor, std::abs(vec_value));
3962 : }
3963 : }
3964 228330 : };
3965 :
3966 : // Compute our scaling factors for the spatial field variables
3967 82245 : for (const auto & elem : _fe_problem.getCurrentAlgebraicElementRange())
3968 310836 : for (const auto i : make_range(system().n_vars()))
3969 229243 : if (_variable_autoscaled[i] && system().variable_type(i).family != SCALAR)
3970 : {
3971 228237 : dof_map.dof_indices(elem, dof_indices, i);
3972 228237 : examine_dof_indices(dof_indices, i);
3973 : }
3974 :
3975 1712 : for (const auto i : make_range(system().n_vars()))
3976 1060 : if (_variable_autoscaled[i] && system().variable_type(i).family == SCALAR)
3977 : {
3978 93 : dof_map.SCALAR_dof_indices(dof_indices, i);
3979 93 : examine_dof_indices(dof_indices, i);
3980 : }
3981 :
3982 652 : if (resid_scaling)
3983 50 : _communicator.max(resid_inverse_scaling_factors);
3984 652 : if (jac_scaling)
3985 602 : _communicator.max(jac_inverse_scaling_factors);
3986 :
3987 652 : if (jac_scaling && resid_scaling)
3988 0 : for (MooseIndex(inverse_scaling_factors) i = 0; i < inverse_scaling_factors.size(); ++i)
3989 : {
3990 : // Be careful not to take log(0)
3991 0 : if (!resid_inverse_scaling_factors[i])
3992 : {
3993 0 : if (!jac_inverse_scaling_factors[i])
3994 0 : inverse_scaling_factors[i] = 1;
3995 : else
3996 0 : inverse_scaling_factors[i] = jac_inverse_scaling_factors[i];
3997 : }
3998 0 : else if (!jac_inverse_scaling_factors[i])
3999 : // We know the resid is not zero
4000 0 : inverse_scaling_factors[i] = resid_inverse_scaling_factors[i];
4001 : else
4002 0 : inverse_scaling_factors[i] =
4003 0 : std::exp(_resid_vs_jac_scaling_param * std::log(resid_inverse_scaling_factors[i]) +
4004 0 : (1 - _resid_vs_jac_scaling_param) * std::log(jac_inverse_scaling_factors[i]));
4005 0 : }
4006 652 : else if (jac_scaling)
4007 602 : inverse_scaling_factors = jac_inverse_scaling_factors;
4008 50 : else if (resid_scaling)
4009 50 : inverse_scaling_factors = resid_inverse_scaling_factors;
4010 : else
4011 0 : mooseError("We shouldn't be calling this routine if we're not performing any scaling");
4012 :
4013 : // We have to make sure that our scaling values are not zero
4014 1688 : for (auto & scaling_factor : inverse_scaling_factors)
4015 1036 : if (scaling_factor == 0)
4016 40 : scaling_factor = 1;
4017 :
4018 : // Now flatten the group scaling factors to the individual variable scaling factors
4019 652 : std::vector<Real> flattened_inverse_scaling_factors(system().n_vars());
4020 1712 : for (const auto i : index_range(flattened_inverse_scaling_factors))
4021 1060 : flattened_inverse_scaling_factors[i] = inverse_scaling_factors[_var_to_group_var[i]];
4022 :
4023 : // Now set the scaling factors for the variables
4024 652 : applyScalingFactors(flattened_inverse_scaling_factors);
4025 652 : if (auto displaced_problem = _fe_problem.getDisplacedProblem().get())
4026 63 : displaced_problem->systemBaseNonlinear(number()).applyScalingFactors(
4027 : flattened_inverse_scaling_factors);
4028 :
4029 652 : _computed_scaling = true;
4030 652 : return true;
4031 652 : }
4032 :
4033 : void
4034 314500 : NonlinearSystemBase::assembleScalingVector()
4035 : {
4036 943500 : if (!hasVector("scaling_factors"))
4037 : // No variables have indicated they need scaling
4038 313950 : return;
4039 :
4040 1100 : auto & scaling_vector = getVector("scaling_factors");
4041 :
4042 550 : const auto & lm_mesh = _mesh.getMesh();
4043 550 : const auto & dof_map = dofMap();
4044 :
4045 550 : const auto & field_variables = _vars[0].fieldVariables();
4046 550 : const auto & scalar_variables = _vars[0].scalars();
4047 :
4048 550 : std::vector<dof_id_type> dof_indices;
4049 :
4050 550 : for (const Elem * const elem :
4051 13672 : as_range(lm_mesh.active_local_elements_begin(), lm_mesh.active_local_elements_end()))
4052 31020 : for (const auto * const field_var : field_variables)
4053 : {
4054 18448 : const auto & factors = field_var->arrayScalingFactor();
4055 36896 : for (const auto i : make_range(field_var->count()))
4056 : {
4057 18448 : dof_map.dof_indices(elem, dof_indices, field_var->number() + i);
4058 76448 : for (const auto dof : dof_indices)
4059 58000 : scaling_vector.set(dof, factors[i]);
4060 : }
4061 550 : }
4062 :
4063 662 : for (const auto * const scalar_var : scalar_variables)
4064 : {
4065 : mooseAssert(scalar_var->count() == 1,
4066 : "Scalar variables should always have only one component.");
4067 112 : dof_map.SCALAR_dof_indices(dof_indices, scalar_var->number());
4068 224 : for (const auto dof : dof_indices)
4069 112 : scaling_vector.set(dof, scalar_var->scalingFactor());
4070 : }
4071 :
4072 : // Parallel assemble
4073 550 : scaling_vector.close();
4074 :
4075 550 : if (auto * displaced_problem = _fe_problem.getDisplacedProblem().get())
4076 : // copy into the corresponding displaced system vector because they should be the exact same
4077 0 : displaced_problem->systemBaseNonlinear(number()).getVector("scaling_factors") = scaling_vector;
4078 550 : }
4079 :
4080 : bool
4081 313848 : NonlinearSystemBase::preSolve()
4082 : {
4083 : // Clear the iteration counters
4084 313848 : _current_l_its.clear();
4085 313848 : _current_nl_its = 0;
4086 :
4087 : // Initialize the solution vector using a predictor and known values from nodal bcs
4088 313848 : setInitialSolution();
4089 :
4090 : // Now that the initial solution has ben set, potentially perform a residual/Jacobian evaluation
4091 : // to determine variable scaling factors
4092 313848 : if (_automatic_scaling)
4093 : {
4094 1112 : const bool scaling_succeeded = computeScaling();
4095 1112 : if (!scaling_succeeded)
4096 0 : return false;
4097 : }
4098 :
4099 : // We do not know a priori what variable a global degree of freedom corresponds to, so we need a
4100 : // map from global dof to scaling factor. We just use a ghosted NumericVector for that mapping
4101 313848 : assembleScalingVector();
4102 :
4103 313848 : return true;
4104 : }
4105 :
4106 : void
4107 5612 : NonlinearSystemBase::destroyColoring()
4108 : {
4109 5612 : if (matrixFromColoring())
4110 117 : LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
4111 5612 : }
4112 :
4113 : FieldSplitPreconditionerBase &
4114 0 : NonlinearSystemBase::getFieldSplitPreconditioner()
4115 : {
4116 0 : if (!_fsp)
4117 0 : mooseError("No field split preconditioner is present for this system");
4118 :
4119 0 : return *_fsp;
4120 : }
|
