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