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 61453 : NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem,
114 : System & sys,
115 61453 : const std::string & name)
116 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
117 61453 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
118 61453 : _sys(sys),
119 61453 : _last_nl_rnorm(0.),
120 61453 : _current_nl_its(0),
121 61453 : _residual_ghosted(NULL),
122 61453 : _Re_time_tag(-1),
123 61453 : _Re_time(NULL),
124 61453 : _Re_non_time_tag(-1),
125 61453 : _Re_non_time(NULL),
126 61453 : _scalar_kernels(/*threaded=*/false),
127 61453 : _nodal_bcs(/*threaded=*/false),
128 61453 : _preset_nodal_bcs(/*threaded=*/false),
129 61453 : _ad_preset_nodal_bcs(/*threaded=*/false),
130 61453 : _splits(/*threaded=*/false),
131 61453 : _increment_vec(NULL),
132 61453 : _use_finite_differenced_preconditioner(false),
133 61453 : _fdcoloring(nullptr),
134 61453 : _have_decomposition(false),
135 61453 : _use_field_split_preconditioner(false),
136 61453 : _add_implicit_geometric_coupling_entries_to_jacobian(false),
137 61453 : _assemble_constraints_separately(false),
138 61453 : _need_residual_ghosted(false),
139 61453 : _debugging_residuals(false),
140 61453 : _doing_dg(false),
141 61453 : _n_iters(0),
142 61453 : _n_linear_iters(0),
143 61453 : _n_residual_evaluations(0),
144 61453 : _final_residual(0.),
145 61453 : _computing_pre_smo_residual(false),
146 61453 : _pre_smo_residual(0),
147 61453 : _initial_residual(0),
148 61453 : _use_pre_smo_residual(false),
149 61453 : _print_all_var_norms(false),
150 61453 : _has_save_in(false),
151 61453 : _has_diag_save_in(false),
152 61453 : _has_nodalbc_save_in(false),
153 61453 : _has_nodalbc_diag_save_in(false),
154 61453 : _computed_scaling(false),
155 61453 : _compute_scaling_once(true),
156 61453 : _resid_vs_jac_scaling_param(0),
157 61453 : _off_diagonals_in_auto_scaling(false),
158 307265 : _auto_scaling_initd(false)
159 : {
160 61453 : getResidualNonTimeVector();
161 : // Don't need to add the matrix - it already exists (for now)
162 61453 : _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 61453 : _fe_problem.addMatrixTag("TIME");
167 :
168 61453 : _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
169 :
170 61453 : _sys.identify_variable_groups(_fe_problem.identifyVariableGroupsInNL());
171 :
172 61453 : if (!_fe_problem.defaultGhosting())
173 : {
174 61371 : auto & dof_map = _sys.get_dof_map();
175 61371 : dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting());
176 61371 : dof_map.set_implicit_neighbor_dofs(false);
177 : }
178 61453 : }
179 :
180 57125 : NonlinearSystemBase::~NonlinearSystemBase() = default;
181 :
182 : void
183 59555 : NonlinearSystemBase::preInit()
184 : {
185 59555 : SolverSystem::preInit();
186 :
187 59555 : if (_fe_problem.hasDampers())
188 159 : setupDampers();
189 :
190 59555 : if (_residual_copy.get())
191 0 : _residual_copy->init(_sys.n_dofs(), false, SERIAL);
192 59555 : }
193 :
194 : void
195 76 : NonlinearSystemBase::turnOffJacobian()
196 : {
197 76 : system().set_basic_system_only();
198 76 : nonlinearSolver()->jacobian = NULL;
199 76 : }
200 :
201 : void
202 58668 : NonlinearSystemBase::initialSetup()
203 : {
204 58668 : TIME_SECTION("nlInitialSetup", 2, "Setting Up Nonlinear System");
205 :
206 58668 : SolverSystem::initialSetup();
207 :
208 : {
209 58668 : TIME_SECTION("kernelsInitialSetup", 2, "Setting Up Kernels/BCs/Constraints");
210 :
211 122837 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
212 : {
213 64169 : _kernels.initialSetup(tid);
214 64169 : _nodal_kernels.initialSetup(tid);
215 64169 : _dirac_kernels.initialSetup(tid);
216 64169 : if (_doing_dg)
217 1198 : _dg_kernels.initialSetup(tid);
218 64169 : _interface_kernels.initialSetup(tid);
219 :
220 64169 : _element_dampers.initialSetup(tid);
221 64169 : _nodal_dampers.initialSetup(tid);
222 64169 : _integrated_bcs.initialSetup(tid);
223 :
224 64169 : if (_fe_problem.haveFV())
225 : {
226 4283 : std::vector<FVElementalKernel *> fv_elemental_kernels;
227 4283 : _fe_problem.theWarehouse()
228 8566 : .query()
229 4283 : .template condition<AttribSystem>("FVElementalKernel")
230 4283 : .template condition<AttribThread>(tid)
231 4283 : .queryInto(fv_elemental_kernels);
232 :
233 8231 : for (auto * fv_kernel : fv_elemental_kernels)
234 3948 : fv_kernel->initialSetup();
235 :
236 4283 : std::vector<FVFluxKernel *> fv_flux_kernels;
237 4283 : _fe_problem.theWarehouse()
238 8566 : .query()
239 4283 : .template condition<AttribSystem>("FVFluxKernel")
240 4283 : .template condition<AttribThread>(tid)
241 4283 : .queryInto(fv_flux_kernels);
242 :
243 9686 : for (auto * fv_kernel : fv_flux_kernels)
244 5403 : fv_kernel->initialSetup();
245 4283 : }
246 : }
247 :
248 58668 : _scalar_kernels.initialSetup();
249 58668 : _constraints.initialSetup();
250 58668 : _general_dampers.initialSetup();
251 58668 : _nodal_bcs.initialSetup();
252 58668 : }
253 :
254 : {
255 58668 : TIME_SECTION("mortarSetup", 2, "Initializing Mortar Interfaces");
256 :
257 123105 : auto create_mortar_functors = [this](const bool displaced)
258 : {
259 : // go over mortar interfaces and construct functors
260 117336 : const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
261 118330 : for (const auto & mortar_interface : mortar_interfaces)
262 : {
263 994 : const auto primary_secondary_boundary_pair = mortar_interface.first;
264 994 : if (!_constraints.hasActiveMortarConstraints(primary_secondary_boundary_pair, displaced))
265 39 : continue;
266 :
267 955 : const auto & mortar_generation_object = mortar_interface.second;
268 :
269 : auto & mortar_constraints =
270 955 : _constraints.getActiveMortarConstraints(primary_secondary_boundary_pair, displaced);
271 :
272 : auto & subproblem = displaced
273 196 : ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
274 1053 : : static_cast<SubProblem &>(_fe_problem);
275 :
276 955 : auto & mortar_functors =
277 : displaced ? _displaced_mortar_functors : _undisplaced_mortar_functors;
278 :
279 955 : mortar_functors.emplace(primary_secondary_boundary_pair,
280 1910 : ComputeMortarFunctor(mortar_constraints,
281 : mortar_generation_object,
282 : subproblem,
283 : _fe_problem,
284 : displaced,
285 955 : subproblem.assembly(0, number())));
286 : }
287 117336 : };
288 :
289 58668 : create_mortar_functors(false);
290 58668 : create_mortar_functors(true);
291 58668 : }
292 :
293 58668 : if (_automatic_scaling)
294 : {
295 445 : if (_off_diagonals_in_auto_scaling)
296 69 : _scaling_matrix = std::make_unique<OffDiagonalScalingMatrix<Number>>(_communicator);
297 : else
298 376 : _scaling_matrix = std::make_unique<DiagonalMatrix<Number>>(_communicator);
299 : }
300 :
301 58668 : if (_preconditioner)
302 14079 : _preconditioner->initialSetup();
303 58668 : }
304 :
305 : void
306 285354 : NonlinearSystemBase::timestepSetup()
307 : {
308 285354 : SolverSystem::timestepSetup();
309 :
310 596148 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
311 : {
312 310794 : _kernels.timestepSetup(tid);
313 310794 : _nodal_kernels.timestepSetup(tid);
314 310794 : _dirac_kernels.timestepSetup(tid);
315 310794 : if (_doing_dg)
316 1775 : _dg_kernels.timestepSetup(tid);
317 310794 : _interface_kernels.timestepSetup(tid);
318 310794 : _element_dampers.timestepSetup(tid);
319 310794 : _nodal_dampers.timestepSetup(tid);
320 310794 : _integrated_bcs.timestepSetup(tid);
321 :
322 310794 : if (_fe_problem.haveFV())
323 : {
324 17387 : std::vector<FVFluxBC *> bcs;
325 17387 : _fe_problem.theWarehouse()
326 34774 : .query()
327 17387 : .template condition<AttribSystem>("FVFluxBC")
328 17387 : .template condition<AttribThread>(tid)
329 17387 : .queryInto(bcs);
330 :
331 17387 : std::vector<FVInterfaceKernel *> iks;
332 17387 : _fe_problem.theWarehouse()
333 34774 : .query()
334 17387 : .template condition<AttribSystem>("FVInterfaceKernel")
335 17387 : .template condition<AttribThread>(tid)
336 17387 : .queryInto(iks);
337 :
338 17387 : std::vector<FVFluxKernel *> kernels;
339 17387 : _fe_problem.theWarehouse()
340 34774 : .query()
341 17387 : .template condition<AttribSystem>("FVFluxKernel")
342 17387 : .template condition<AttribThread>(tid)
343 17387 : .queryInto(kernels);
344 :
345 34378 : for (auto * bc : bcs)
346 16991 : bc->timestepSetup();
347 17785 : for (auto * ik : iks)
348 398 : ik->timestepSetup();
349 38501 : for (auto * kernel : kernels)
350 21114 : kernel->timestepSetup();
351 17387 : }
352 : }
353 285354 : _scalar_kernels.timestepSetup();
354 285354 : _constraints.timestepSetup();
355 285354 : _general_dampers.timestepSetup();
356 285354 : _nodal_bcs.timestepSetup();
357 285354 : }
358 :
359 : void
360 1911830 : NonlinearSystemBase::customSetup(const ExecFlagType & exec_type)
361 : {
362 1911830 : SolverSystem::customSetup(exec_type);
363 :
364 3995534 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
365 : {
366 2083704 : _kernels.customSetup(exec_type, tid);
367 2083704 : _nodal_kernels.customSetup(exec_type, tid);
368 2083704 : _dirac_kernels.customSetup(exec_type, tid);
369 2083704 : if (_doing_dg)
370 12147 : _dg_kernels.customSetup(exec_type, tid);
371 2083704 : _interface_kernels.customSetup(exec_type, tid);
372 2083704 : _element_dampers.customSetup(exec_type, tid);
373 2083704 : _nodal_dampers.customSetup(exec_type, tid);
374 2083704 : _integrated_bcs.customSetup(exec_type, tid);
375 :
376 2083704 : if (_fe_problem.haveFV())
377 : {
378 127728 : std::vector<FVFluxBC *> bcs;
379 127728 : _fe_problem.theWarehouse()
380 255456 : .query()
381 127728 : .template condition<AttribSystem>("FVFluxBC")
382 127728 : .template condition<AttribThread>(tid)
383 127728 : .queryInto(bcs);
384 :
385 127728 : std::vector<FVInterfaceKernel *> iks;
386 127728 : _fe_problem.theWarehouse()
387 255456 : .query()
388 127728 : .template condition<AttribSystem>("FVInterfaceKernel")
389 127728 : .template condition<AttribThread>(tid)
390 127728 : .queryInto(iks);
391 :
392 127728 : std::vector<FVFluxKernel *> kernels;
393 127728 : _fe_problem.theWarehouse()
394 255456 : .query()
395 127728 : .template condition<AttribSystem>("FVFluxKernel")
396 127728 : .template condition<AttribThread>(tid)
397 127728 : .queryInto(kernels);
398 :
399 230660 : for (auto * bc : bcs)
400 102932 : bc->customSetup(exec_type);
401 171947 : for (auto * ik : iks)
402 44219 : ik->customSetup(exec_type);
403 304236 : for (auto * kernel : kernels)
404 176508 : kernel->customSetup(exec_type);
405 127728 : }
406 : }
407 1911830 : _scalar_kernels.customSetup(exec_type);
408 1911830 : _constraints.customSetup(exec_type);
409 1911830 : _general_dampers.customSetup(exec_type);
410 1911830 : _nodal_bcs.customSetup(exec_type);
411 1911830 : }
412 :
413 : void
414 341883 : NonlinearSystemBase::setupDM()
415 : {
416 341883 : if (haveFieldSplitPreconditioner())
417 105 : Moose::PetscSupport::petscSetupDM(*this, _decomposition_split);
418 341883 : }
419 :
420 : void
421 60961 : 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 60961 : if (splits.size() && splits.size() != 1)
425 0 : mooseError("Only a single top-level split is allowed in a Problem's decomposition.");
426 :
427 60961 : if (splits.size())
428 : {
429 113 : _decomposition_split = splits[0];
430 113 : _have_decomposition = true;
431 : }
432 : else
433 60848 : _have_decomposition = false;
434 60961 : }
435 :
436 : void
437 113 : NonlinearSystemBase::setupFieldDecomposition()
438 : {
439 113 : if (!_have_decomposition)
440 0 : return;
441 :
442 113 : std::shared_ptr<Split> top_split = getSplit(_decomposition_split);
443 113 : top_split->setup(*this);
444 113 : }
445 :
446 : void
447 81806 : NonlinearSystemBase::addKernel(const std::string & kernel_name,
448 : const std::string & name,
449 : InputParameters & parameters)
450 : {
451 171026 : 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 89420 : _factory.create<KernelBase>(kernel_name, name, parameters, tid);
456 89244 : _kernels.addObject(kernel, tid);
457 89224 : postAddResidualObject(*kernel);
458 : // Add to theWarehouse, a centralized storage for all moose objects
459 89220 : _fe_problem.theWarehouse().add(kernel);
460 89220 : }
461 :
462 81606 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
463 130 : _has_save_in = true;
464 81606 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
465 104 : _has_diag_save_in = true;
466 81606 : }
467 :
468 : void
469 769 : NonlinearSystemBase::addHDGKernel(const std::string & kernel_name,
470 : const std::string & name,
471 : InputParameters & parameters)
472 : {
473 1543 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
474 : {
475 : // Create the kernel object via the factory and add to warehouse
476 774 : auto kernel = _factory.create<HDGKernel>(kernel_name, name, parameters, tid);
477 774 : _kernels.addObject(kernel, tid);
478 774 : _hybridized_kernels.addObject(kernel, tid);
479 : // Add to theWarehouse, a centralized storage for all moose objects
480 774 : _fe_problem.theWarehouse().add(kernel);
481 774 : postAddResidualObject(*kernel);
482 774 : }
483 769 : }
484 :
485 : void
486 570 : NonlinearSystemBase::addNodalKernel(const std::string & kernel_name,
487 : const std::string & name,
488 : InputParameters & parameters)
489 : {
490 1241 : 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 671 : _factory.create<NodalKernelBase>(kernel_name, name, parameters, tid);
495 671 : _nodal_kernels.addObject(kernel, tid);
496 : // Add to theWarehouse, a centralized storage for all moose objects
497 671 : _fe_problem.theWarehouse().add(kernel);
498 671 : postAddResidualObject(*kernel);
499 671 : }
500 :
501 570 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
502 0 : _has_save_in = true;
503 570 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
504 0 : _has_diag_save_in = true;
505 570 : }
506 :
507 : void
508 1525 : NonlinearSystemBase::addScalarKernel(const std::string & kernel_name,
509 : const std::string & name,
510 : InputParameters & parameters)
511 : {
512 : std::shared_ptr<ScalarKernelBase> kernel =
513 1525 : _factory.create<ScalarKernelBase>(kernel_name, name, parameters);
514 1517 : postAddResidualObject(*kernel);
515 : // Add to theWarehouse, a centralized storage for all moose objects
516 1517 : _fe_problem.theWarehouse().add(kernel);
517 1517 : _scalar_kernels.addObject(kernel);
518 1517 : }
519 :
520 : void
521 81180 : NonlinearSystemBase::addBoundaryCondition(const std::string & bc_name,
522 : const std::string & name,
523 : InputParameters & parameters)
524 : {
525 : // ThreadID
526 81180 : THREAD_ID tid = 0;
527 :
528 : // Create the object
529 : std::shared_ptr<BoundaryCondition> bc =
530 81180 : _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
531 81127 : postAddResidualObject(*bc);
532 :
533 : // Active BoundaryIDs for the object
534 81127 : const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
535 81127 : auto bc_var = dynamic_cast<const MooseVariableFieldBase *>(&bc->variable());
536 81127 : _vars[tid].addBoundaryVar(boundary_ids, bc_var);
537 :
538 : // Cast to the various types of BCs
539 81127 : std::shared_ptr<NodalBCBase> nbc = std::dynamic_pointer_cast<NodalBCBase>(bc);
540 81127 : std::shared_ptr<IntegratedBCBase> ibc = std::dynamic_pointer_cast<IntegratedBCBase>(bc);
541 :
542 : // NodalBCBase
543 81127 : if (nbc)
544 : {
545 72300 : 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 72296 : _nodal_bcs.addObject(nbc);
553 : // Add to theWarehouse, a centralized storage for all moose objects
554 72296 : _fe_problem.theWarehouse().add(nbc);
555 72296 : _vars[tid].addBoundaryVars(boundary_ids, nbc->getCoupledVars());
556 :
557 72296 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
558 99 : _has_nodalbc_save_in = true;
559 72296 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
560 21 : _has_nodalbc_diag_save_in = true;
561 :
562 : // DirichletBCs that are preset
563 72296 : std::shared_ptr<DirichletBCBase> dbc = std::dynamic_pointer_cast<DirichletBCBase>(bc);
564 72296 : if (dbc && dbc->preset())
565 62123 : _preset_nodal_bcs.addObject(dbc);
566 :
567 72296 : std::shared_ptr<ADDirichletBCBase> addbc = std::dynamic_pointer_cast<ADDirichletBCBase>(bc);
568 72296 : if (addbc && addbc->preset())
569 1636 : _ad_preset_nodal_bcs.addObject(addbc);
570 72296 : }
571 :
572 : // IntegratedBCBase
573 8827 : else if (ibc)
574 : {
575 8827 : _integrated_bcs.addObject(ibc, tid);
576 : // Add to theWarehouse, a centralized storage for all moose objects
577 8827 : _fe_problem.theWarehouse().add(ibc);
578 8827 : _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
579 :
580 8827 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
581 86 : _has_save_in = true;
582 8827 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
583 60 : _has_diag_save_in = true;
584 :
585 9487 : for (tid = 1; tid < libMesh::n_threads(); tid++)
586 : {
587 : // Create the object
588 660 : bc = _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
589 :
590 : // Give users opportunity to set some parameters
591 660 : postAddResidualObject(*bc);
592 :
593 : // Active BoundaryIDs for the object
594 660 : const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
595 660 : _vars[tid].addBoundaryVar(boundary_ids, bc_var);
596 :
597 660 : ibc = std::static_pointer_cast<IntegratedBCBase>(bc);
598 :
599 660 : _integrated_bcs.addObject(ibc, tid);
600 660 : _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
601 : }
602 : }
603 :
604 : else
605 0 : mooseError("Unknown BoundaryCondition type for object named ", bc->name());
606 81123 : }
607 :
608 : void
609 1701 : NonlinearSystemBase::addConstraint(const std::string & c_name,
610 : const std::string & name,
611 : InputParameters & parameters)
612 : {
613 1701 : std::shared_ptr<Constraint> constraint = _factory.create<Constraint>(c_name, name, parameters);
614 1693 : _constraints.addObject(constraint);
615 1693 : postAddResidualObject(*constraint);
616 :
617 1693 : if (constraint && constraint->addCouplingEntriesToJacobian())
618 1680 : addImplicitGeometricCouplingEntriesToJacobian(true);
619 1693 : }
620 :
621 : void
622 915 : NonlinearSystemBase::addDiracKernel(const std::string & kernel_name,
623 : const std::string & name,
624 : InputParameters & parameters)
625 : {
626 1906 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
627 : {
628 : std::shared_ptr<DiracKernelBase> kernel =
629 999 : _factory.create<DiracKernelBase>(kernel_name, name, parameters, tid);
630 991 : postAddResidualObject(*kernel);
631 991 : _dirac_kernels.addObject(kernel, tid);
632 : // Add to theWarehouse, a centralized storage for all moose objects
633 991 : _fe_problem.theWarehouse().add(kernel);
634 991 : }
635 907 : }
636 :
637 : void
638 1402 : NonlinearSystemBase::addDGKernel(std::string dg_kernel_name,
639 : const std::string & name,
640 : InputParameters & parameters)
641 : {
642 2923 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
643 : {
644 1521 : auto dg_kernel = _factory.create<DGKernelBase>(dg_kernel_name, name, parameters, tid);
645 1521 : _dg_kernels.addObject(dg_kernel, tid);
646 : // Add to theWarehouse, a centralized storage for all moose objects
647 1521 : _fe_problem.theWarehouse().add(dg_kernel);
648 1521 : postAddResidualObject(*dg_kernel);
649 1521 : }
650 :
651 1402 : _doing_dg = true;
652 :
653 1402 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
654 52 : _has_save_in = true;
655 1402 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
656 39 : _has_diag_save_in = true;
657 1402 : }
658 :
659 : void
660 907 : NonlinearSystemBase::addInterfaceKernel(std::string interface_kernel_name,
661 : const std::string & name,
662 : InputParameters & parameters)
663 : {
664 1889 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
665 : {
666 : std::shared_ptr<InterfaceKernelBase> interface_kernel =
667 982 : _factory.create<InterfaceKernelBase>(interface_kernel_name, name, parameters, tid);
668 982 : postAddResidualObject(*interface_kernel);
669 :
670 982 : const std::set<BoundaryID> & boundary_ids = interface_kernel->boundaryIDs();
671 982 : auto ik_var = dynamic_cast<const MooseVariableFieldBase *>(&interface_kernel->variable());
672 982 : _vars[tid].addBoundaryVar(boundary_ids, ik_var);
673 :
674 982 : _interface_kernels.addObject(interface_kernel, tid);
675 : // Add to theWarehouse, a centralized storage for all moose objects
676 982 : _fe_problem.theWarehouse().add(interface_kernel);
677 982 : _vars[tid].addBoundaryVars(boundary_ids, interface_kernel->getCoupledVars());
678 982 : }
679 907 : }
680 :
681 : void
682 177 : NonlinearSystemBase::addDamper(const std::string & damper_name,
683 : const std::string & name,
684 : InputParameters & parameters)
685 : {
686 328 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
687 : {
688 218 : std::shared_ptr<Damper> damper = _factory.create<Damper>(damper_name, name, parameters, tid);
689 :
690 : // Attempt to cast to the damper types
691 218 : std::shared_ptr<ElementDamper> ed = std::dynamic_pointer_cast<ElementDamper>(damper);
692 218 : std::shared_ptr<NodalDamper> nd = std::dynamic_pointer_cast<NodalDamper>(damper);
693 218 : std::shared_ptr<GeneralDamper> gd = std::dynamic_pointer_cast<GeneralDamper>(damper);
694 :
695 218 : if (gd)
696 : {
697 67 : _general_dampers.addObject(gd);
698 67 : 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 419 : }
707 177 : }
708 :
709 : void
710 407 : NonlinearSystemBase::addSplit(const std::string & split_name,
711 : const std::string & name,
712 : InputParameters & parameters)
713 : {
714 407 : std::shared_ptr<Split> split = _factory.create<Split>(split_name, name, parameters);
715 407 : _splits.addObject(split);
716 : // Add to theWarehouse, a centralized storage for all moose objects
717 407 : _fe_problem.theWarehouse().add(split);
718 407 : }
719 :
720 : std::shared_ptr<Split>
721 407 : NonlinearSystemBase::getSplit(const std::string & name)
722 : {
723 407 : return _splits.getActiveObject(name);
724 : }
725 :
726 : bool
727 311354 : NonlinearSystemBase::shouldEvaluatePreSMOResidual() const
728 : {
729 311354 : if (_fe_problem.solverParams(number())._type == Moose::ST_LINEAR)
730 10637 : 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 300717 : if (_app.parameters().get<bool>("use_legacy_initial_residual_evaluation_behavior"))
739 0 : return true;
740 :
741 300717 : return _use_pre_smo_residual;
742 : }
743 :
744 : Real
745 385248 : NonlinearSystemBase::referenceResidual() const
746 : {
747 385248 : return usePreSMOResidual() ? preSMOResidual() : initialResidual();
748 : }
749 :
750 : Real
751 171 : NonlinearSystemBase::preSMOResidual() const
752 : {
753 171 : if (!shouldEvaluatePreSMOResidual())
754 0 : mooseError("pre-SMO residual is requested but not evaluated.");
755 :
756 171 : return _pre_smo_residual;
757 : }
758 :
759 : Real
760 385360 : NonlinearSystemBase::initialResidual() const
761 : {
762 385360 : return _initial_residual;
763 : }
764 :
765 : void
766 308757 : NonlinearSystemBase::setInitialResidual(Real r)
767 : {
768 308757 : _initial_residual = r;
769 308757 : }
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 131 : NonlinearSystemBase::computeResidualTag(NumericVector<Number> & residual, TagID tag_id)
783 : {
784 131 : _nl_vector_tags.clear();
785 131 : _nl_vector_tags.insert(tag_id);
786 131 : _nl_vector_tags.insert(residualVectorTag());
787 :
788 131 : associateVectorToTag(residual, residualVectorTag());
789 :
790 131 : computeResidualTags(_nl_vector_tags);
791 :
792 131 : disassociateVectorFromTag(residual, residualVectorTag());
793 131 : }
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 3278015 : NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
805 : {
806 : parallel_object_only();
807 :
808 3278015 : TIME_SECTION("nl::computeResidualTags", 5);
809 :
810 3278015 : _fe_problem.setCurrentNonlinearSystem(number());
811 3278015 : _fe_problem.setCurrentlyComputingResidual(true);
812 :
813 3278015 : bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
814 :
815 3278015 : _n_residual_evaluations++;
816 :
817 : // not suppose to do anythin on matrix
818 3278015 : deactiveAllMatrixTags();
819 :
820 3278015 : FloatingPointExceptionGuard fpe_guard(_app);
821 :
822 3278015 : 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 3278015 : zeroTaggedVectors(tags);
833 3278015 : computeResidualInternal(tags);
834 3277853 : closeTaggedVectors(tags);
835 :
836 3277853 : if (required_residual)
837 : {
838 3245971 : auto & residual = getVector(residualVectorTag());
839 3245971 : if (!_time_integrators.empty())
840 : {
841 5709381 : for (auto & ti : _time_integrators)
842 2863116 : ti->postResidual(residual);
843 : }
844 : else
845 399706 : residual += *_Re_non_time;
846 3245971 : residual.close();
847 : }
848 3277853 : if (_fe_problem.computingScalingResidual())
849 : // We don't want to do nodal bcs or anything else
850 50 : return;
851 :
852 3277803 : computeNodalBCs(tags);
853 3277803 : closeTaggedVectors(tags);
854 :
855 : // If we are debugging residuals we need one more assignment to have the ghosted copy up to
856 : // date
857 3277803 : 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 3277803 : if (_has_nodalbc_save_in)
866 171 : _fe_problem.getAuxiliarySystem().solution().close();
867 3277803 : if (hasSaveIn())
868 309 : _fe_problem.getAuxiliarySystem().update();
869 : }
870 121 : 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 121 : }
876 :
877 : // not supposed to do anything on matrix
878 3277924 : activeAllMatrixTags();
879 :
880 3277924 : _fe_problem.setCurrentlyComputingResidual(false);
881 3278024 : }
882 :
883 : void
884 3467 : NonlinearSystemBase::computeResidualAndJacobianTags(const std::set<TagID> & vector_tags,
885 : const std::set<TagID> & matrix_tags)
886 : {
887 : const bool required_residual =
888 3467 : vector_tags.find(residualVectorTag()) == vector_tags.end() ? false : true;
889 :
890 : try
891 : {
892 3467 : zeroTaggedVectors(vector_tags);
893 3467 : computeResidualAndJacobianInternal(vector_tags, matrix_tags);
894 3467 : closeTaggedVectors(vector_tags);
895 3467 : closeTaggedMatrices(matrix_tags);
896 :
897 3467 : if (required_residual)
898 : {
899 3467 : auto & residual = getVector(residualVectorTag());
900 3467 : if (!_time_integrators.empty())
901 : {
902 3818 : for (auto & ti : _time_integrators)
903 1909 : ti->postResidual(residual);
904 : }
905 : else
906 1558 : residual += *_Re_non_time;
907 3467 : residual.close();
908 : }
909 :
910 3467 : computeNodalBCsResidualAndJacobian();
911 3467 : closeTaggedVectors(vector_tags);
912 3467 : 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 3467 : }
921 :
922 : void
923 256795 : NonlinearSystemBase::onTimestepBegin()
924 : {
925 514308 : for (auto & ti : _time_integrators)
926 257513 : ti->preSolve();
927 256795 : if (_predictor.get())
928 387 : _predictor->timestepSetup();
929 256795 : }
930 :
931 : void
932 312489 : NonlinearSystemBase::setInitialSolution()
933 : {
934 312489 : deactiveAllMatrixTags();
935 :
936 312489 : NumericVector<Number> & initial_solution(solution());
937 312489 : if (_predictor.get())
938 : {
939 387 : if (_predictor->shouldApply())
940 : {
941 223 : TIME_SECTION("applyPredictor", 2, "Applying Predictor");
942 :
943 223 : _predictor->apply(initial_solution);
944 223 : _fe_problem.predictorCleanup(initial_solution);
945 223 : }
946 : else
947 164 : _console << " Skipping predictor this step" << std::endl;
948 : }
949 :
950 : // do nodal BC
951 : {
952 312489 : TIME_SECTION("initialBCs", 2, "Applying BCs To Initial Condition");
953 :
954 312489 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
955 15343141 : for (const auto & bnode : bnd_nodes)
956 : {
957 15030652 : BoundaryID boundary_id = bnode->_bnd_id;
958 15030652 : Node * node = bnode->_node;
959 :
960 15030652 : if (node->processor_id() == processor_id())
961 : {
962 : // reinit variables in nodes
963 11521200 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
964 :
965 11521200 : if (_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
966 : {
967 3621993 : const auto & preset_bcs = _preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
968 7892390 : for (const auto & preset_bc : preset_bcs)
969 4270397 : preset_bc->computeValue(initial_solution);
970 : }
971 11521200 : if (_ad_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
972 : {
973 36720 : const auto & preset_bcs_res = _ad_preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
974 77278 : for (const auto & preset_bc : preset_bcs_res)
975 40558 : preset_bc->computeValue(initial_solution);
976 : }
977 : }
978 : }
979 312489 : }
980 :
981 312489 : _sys.solution->close();
982 312489 : update();
983 :
984 : // Set constraint secondary values
985 312489 : setConstraintSecondaryValues(initial_solution, false);
986 :
987 312489 : if (_fe_problem.getDisplacedProblem())
988 35724 : setConstraintSecondaryValues(initial_solution, true);
989 312489 : }
990 :
991 : void
992 121 : NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
993 : {
994 121 : _predictor = predictor;
995 121 : }
996 :
997 : void
998 5983481 : NonlinearSystemBase::subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
999 : {
1000 5983481 : SolverSystem::subdomainSetup();
1001 :
1002 5983481 : _kernels.subdomainSetup(subdomain, tid);
1003 5983481 : _nodal_kernels.subdomainSetup(subdomain, tid);
1004 5983481 : _element_dampers.subdomainSetup(subdomain, tid);
1005 5983481 : _nodal_dampers.subdomainSetup(subdomain, tid);
1006 5983481 : }
1007 :
1008 : NumericVector<Number> &
1009 30734 : NonlinearSystemBase::getResidualTimeVector()
1010 : {
1011 30734 : if (!_Re_time)
1012 : {
1013 30628 : _Re_time_tag = _fe_problem.addVectorTag("TIME");
1014 :
1015 : // Most applications don't need the expense of ghosting
1016 30628 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1017 30628 : _Re_time = &addVector(_Re_time_tag, false, ptype);
1018 : }
1019 106 : 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 30734 : return *_Re_time;
1029 : }
1030 :
1031 : NumericVector<Number> &
1032 92394 : NonlinearSystemBase::getResidualNonTimeVector()
1033 : {
1034 92394 : if (!_Re_non_time)
1035 : {
1036 61453 : _Re_non_time_tag = _fe_problem.addVectorTag("NONTIME");
1037 :
1038 : // Most applications don't need the expense of ghosting
1039 61453 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1040 61453 : _Re_non_time = &addVector(_Re_non_time_tag, false, ptype);
1041 : }
1042 30941 : 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 92394 : 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 23285 : NonlinearSystemBase::enforceNodalConstraintsResidual(NumericVector<Number> & residual)
1073 : {
1074 23285 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1075 23285 : residual.close();
1076 23285 : if (_constraints.hasActiveNodalConstraints())
1077 : {
1078 602 : const auto & ncs = _constraints.getActiveNodalConstraints();
1079 1204 : for (const auto & nc : ncs)
1080 : {
1081 602 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1082 602 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1083 :
1084 602 : if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1085 : {
1086 602 : _fe_problem.reinitNodes(primary_node_ids, tid);
1087 602 : _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1088 602 : nc->computeResidual(residual);
1089 : }
1090 : }
1091 602 : _fe_problem.addCachedResidualDirectly(residual, tid);
1092 602 : residual.close();
1093 : }
1094 23285 : }
1095 :
1096 : void
1097 3617 : NonlinearSystemBase::enforceNodalConstraintsJacobian()
1098 : {
1099 3617 : if (!hasMatrix(systemMatrixTag()))
1100 0 : mooseError(" A system matrix is required");
1101 :
1102 3617 : auto & jacobian = getMatrix(systemMatrixTag());
1103 3617 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1104 :
1105 3617 : if (_constraints.hasActiveNodalConstraints())
1106 : {
1107 121 : const auto & ncs = _constraints.getActiveNodalConstraints();
1108 242 : for (const auto & nc : ncs)
1109 : {
1110 121 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1111 121 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1112 :
1113 121 : if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1114 : {
1115 121 : _fe_problem.reinitNodes(primary_node_ids, tid);
1116 121 : _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1117 121 : nc->computeJacobian(jacobian);
1118 : }
1119 : }
1120 121 : _fe_problem.addCachedJacobian(tid);
1121 : }
1122 3617 : }
1123 :
1124 : void
1125 108178 : NonlinearSystemBase::reinitNodeFace(const Node & secondary_node,
1126 : const BoundaryID secondary_boundary,
1127 : const PenetrationInfo & info,
1128 : const bool displaced)
1129 : {
1130 115408 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1131 165882 : : static_cast<SubProblem &>(_fe_problem);
1132 :
1133 108178 : const Elem * primary_elem = info._elem;
1134 108178 : unsigned int primary_side = info._side_num;
1135 108178 : std::vector<Point> points;
1136 108178 : 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 108178 : _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 108178 : _fe_problem.prepareAssembly(0);
1150 :
1151 108178 : _fe_problem.setNeighborSubdomainID(primary_elem, 0);
1152 :
1153 : //
1154 : // Reinit material on undisplaced mesh
1155 : //
1156 :
1157 : const Elem * const undisplaced_primary_elem =
1158 108178 : 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 108178 : if (displaced)
1163 : {
1164 : const Point reference_point =
1165 57704 : FEMap::inverse_map(primary_elem->dim(), primary_elem, points[0]);
1166 57704 : 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 50474 : return points[0];
1172 108178 : }();
1173 :
1174 108178 : _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 108178 : _fe_problem.reinitMaterialsNeighbor(primary_elem->subdomain_id(), 0, /*swap_stateful=*/false);
1180 :
1181 : // Reinit points for constraint enforcement
1182 108178 : if (displaced)
1183 57704 : subproblem.reinitNeighborPhys(primary_elem, primary_side, points, 0);
1184 108178 : }
1185 :
1186 : void
1187 348213 : 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 71448 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1195 383937 : : static_cast<SubProblem &>(_fe_problem);
1196 348213 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1197 :
1198 348213 : bool constraints_applied = false;
1199 :
1200 408434 : for (const auto & it : penetration_locators)
1201 : {
1202 60221 : PenetrationLocator & pen_loc = *(it.second);
1203 :
1204 60221 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1205 :
1206 60221 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1207 60221 : BoundaryID primary_boundary = pen_loc._primary_boundary;
1208 :
1209 60221 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1210 : {
1211 : const auto & constraints =
1212 2977 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1213 2977 : std::unordered_set<unsigned int> needed_mat_props;
1214 5954 : for (const auto & constraint : constraints)
1215 : {
1216 2977 : const auto & mp_deps = constraint->getMatPropDependencies();
1217 2977 : needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
1218 : }
1219 2977 : _fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/0);
1220 :
1221 42140 : for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1222 : {
1223 39163 : dof_id_type secondary_node_num = secondary_nodes[i];
1224 39163 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1225 :
1226 39163 : if (secondary_node.processor_id() == processor_id())
1227 : {
1228 30944 : if (pen_loc._penetration_info[secondary_node_num])
1229 : {
1230 30944 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1231 :
1232 30944 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1233 :
1234 61888 : for (const auto & nfc : constraints)
1235 : {
1236 30944 : if (nfc->isExplicitConstraint())
1237 28816 : 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 4256 : if (nfc->secondaryBoundary() != secondary_boundary ||
1243 2128 : nfc->primaryBoundary() != primary_boundary)
1244 0 : continue;
1245 :
1246 2128 : if (nfc->shouldApply())
1247 : {
1248 2128 : constraints_applied = true;
1249 2128 : nfc->computeSecondaryValue(solution);
1250 : }
1251 :
1252 2128 : 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 2977 : }
1266 : }
1267 :
1268 : // go over NodeELemConstraints
1269 348213 : std::set<dof_id_type> unique_secondary_node_ids;
1270 :
1271 781461 : for (const auto & secondary_id : _mesh.meshSubdomains())
1272 : {
1273 1155732 : for (const auto & primary_id : _mesh.meshSubdomains())
1274 : {
1275 722484 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1276 : {
1277 : const auto & constraints =
1278 180 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1279 :
1280 : // get unique set of ids of all nodes on current block
1281 180 : unique_secondary_node_ids.clear();
1282 180 : const MeshBase & meshhelper = _mesh.getMesh();
1283 360 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1284 17340 : meshhelper.active_subdomain_elements_end(secondary_id)))
1285 : {
1286 49140 : for (auto & n : elem->node_ref_range())
1287 40740 : unique_secondary_node_ids.insert(n.id());
1288 180 : }
1289 :
1290 14040 : for (auto secondary_node_id : unique_secondary_node_ids)
1291 : {
1292 13860 : Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1293 :
1294 : // check if secondary node is on current processor
1295 13860 : if (secondary_node.processor_id() == processor_id())
1296 : {
1297 : // This reinits the variables that exist on the secondary node
1298 11088 : _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 11088 : _fe_problem.prepareAssembly(0);
1303 :
1304 23584 : for (const auto & nec : constraints)
1305 : {
1306 12496 : if (nec->shouldApply())
1307 : {
1308 5680 : constraints_applied = true;
1309 5680 : nec->computeSecondaryValue(solution);
1310 : }
1311 : }
1312 : }
1313 : }
1314 : }
1315 : }
1316 : }
1317 :
1318 : // See if constraints were applied anywhere
1319 348213 : _communicator.max(constraints_applied);
1320 :
1321 348213 : if (constraints_applied)
1322 : {
1323 719 : solution.close();
1324 719 : update();
1325 : }
1326 348213 : }
1327 :
1328 : void
1329 28101 : NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool displaced)
1330 : {
1331 : // Make sure the residual is in a good state
1332 28101 : 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 9632 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1339 32917 : : static_cast<SubProblem &>(_fe_problem);
1340 28101 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1341 :
1342 : bool constraints_applied;
1343 28101 : bool residual_has_inserted_values = false;
1344 28101 : if (!_assemble_constraints_separately)
1345 28101 : constraints_applied = false;
1346 52167 : for (const auto & it : penetration_locators)
1347 : {
1348 24066 : 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 24066 : PenetrationLocator & pen_loc = *(it.second);
1355 :
1356 24066 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1357 :
1358 24066 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1359 24066 : BoundaryID primary_boundary = pen_loc._primary_boundary;
1360 :
1361 24066 : bool has_writable_variables(false);
1362 :
1363 24066 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1364 : {
1365 : const auto & constraints =
1366 14004 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1367 :
1368 99926 : for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1369 : {
1370 85922 : dof_id_type secondary_node_num = secondary_nodes[i];
1371 85922 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1372 :
1373 85922 : if (secondary_node.processor_id() == processor_id())
1374 : {
1375 73050 : if (pen_loc._penetration_info[secondary_node_num])
1376 : {
1377 73050 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1378 :
1379 73050 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1380 :
1381 146100 : 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 146100 : if (nfc->secondaryBoundary() != secondary_boundary ||
1388 73050 : nfc->primaryBoundary() != primary_boundary)
1389 0 : continue;
1390 :
1391 73050 : if (nfc->shouldApply())
1392 : {
1393 44234 : constraints_applied = true;
1394 44234 : nfc->computeResidual();
1395 :
1396 44234 : 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 44234 : const auto & secondary_var = nfc->variable();
1403 44234 : 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 44234 : std::vector<Number> values = {nfc->secondaryResidual()};
1414 44234 : residual.insert(values, secondary_dofs);
1415 44234 : residual_has_inserted_values = true;
1416 44234 : }
1417 : else
1418 0 : _fe_problem.cacheResidual(0);
1419 44234 : _fe_problem.cacheResidualNeighbor(0);
1420 : }
1421 73050 : if (nfc->hasWritableCoupledVariables())
1422 : {
1423 28816 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1424 28816 : has_writable_variables = true;
1425 57632 : for (auto * var : nfc->getWritableCoupledVariables())
1426 : {
1427 28816 : if (var->isNodalDefined())
1428 28816 : var->insert(_fe_problem.getAuxiliarySystem().solution());
1429 : }
1430 28816 : }
1431 : }
1432 : }
1433 : }
1434 : }
1435 : }
1436 24066 : _communicator.max(has_writable_variables);
1437 :
1438 24066 : 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 2401 : _fe_problem.getAuxiliarySystem().solution().close();
1444 2401 : _fe_problem.getAuxiliarySystem().system().update();
1445 2401 : solutionOld().close();
1446 : }
1447 :
1448 24066 : 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 28101 : if (!_assemble_constraints_separately)
1475 : {
1476 28101 : _communicator.max(constraints_applied);
1477 :
1478 28101 : 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 10864 : _communicator.max(residual_has_inserted_values);
1483 10864 : if (residual_has_inserted_values)
1484 10864 : residual.close();
1485 :
1486 10864 : _fe_problem.addCachedResidualDirectly(residual, 0);
1487 10864 : residual.close();
1488 :
1489 10864 : if (_need_residual_ghosted)
1490 10864 : *_residual_ghosted = residual;
1491 : }
1492 : }
1493 :
1494 : // go over element-element constraint interface
1495 28101 : THREAD_ID tid = 0;
1496 28101 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
1497 28101 : 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 28101 : std::set<dof_id_type> unique_secondary_node_ids;
1543 :
1544 28101 : constraints_applied = false;
1545 28101 : residual_has_inserted_values = false;
1546 28101 : bool has_writable_variables = false;
1547 96096 : for (const auto & secondary_id : _mesh.meshSubdomains())
1548 : {
1549 288326 : for (const auto & primary_id : _mesh.meshSubdomains())
1550 : {
1551 220331 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1552 : {
1553 : const auto & constraints =
1554 360 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1555 :
1556 : // get unique set of ids of all nodes on current block
1557 360 : unique_secondary_node_ids.clear();
1558 360 : const MeshBase & meshhelper = _mesh.getMesh();
1559 720 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1560 34680 : meshhelper.active_subdomain_elements_end(secondary_id)))
1561 : {
1562 98280 : for (auto & n : elem->node_ref_range())
1563 81480 : unique_secondary_node_ids.insert(n.id());
1564 360 : }
1565 :
1566 28080 : for (auto secondary_node_id : unique_secondary_node_ids)
1567 : {
1568 27720 : Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1569 : // check if secondary node is on current processor
1570 27720 : if (secondary_node.processor_id() == processor_id())
1571 : {
1572 : // This reinits the variables that exist on the secondary node
1573 22176 : _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 22176 : _fe_problem.prepareAssembly(0);
1578 :
1579 47168 : for (const auto & nec : constraints)
1580 : {
1581 24992 : if (nec->shouldApply())
1582 : {
1583 11360 : constraints_applied = true;
1584 11360 : nec->computeResidual();
1585 :
1586 11360 : if (nec->overwriteSecondaryResidual())
1587 : {
1588 0 : _fe_problem.setResidual(residual, 0);
1589 0 : residual_has_inserted_values = true;
1590 : }
1591 : else
1592 11360 : _fe_problem.cacheResidual(0);
1593 11360 : _fe_problem.cacheResidualNeighbor(0);
1594 : }
1595 24992 : if (nec->hasWritableCoupledVariables())
1596 : {
1597 352 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1598 352 : has_writable_variables = true;
1599 704 : for (auto * var : nec->getWritableCoupledVariables())
1600 : {
1601 352 : if (var->isNodalDefined())
1602 352 : var->insert(_fe_problem.getAuxiliarySystem().solution());
1603 : }
1604 352 : }
1605 : }
1606 22176 : _fe_problem.addCachedResidual(0);
1607 : }
1608 : }
1609 : }
1610 : }
1611 : }
1612 28101 : _communicator.max(constraints_applied);
1613 :
1614 28101 : 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 360 : _communicator.max(residual_has_inserted_values);
1619 360 : if (residual_has_inserted_values)
1620 0 : residual.close();
1621 :
1622 360 : _fe_problem.addCachedResidualDirectly(residual, 0);
1623 360 : residual.close();
1624 :
1625 360 : if (_need_residual_ghosted)
1626 360 : *_residual_ghosted = residual;
1627 : }
1628 28101 : _communicator.max(has_writable_variables);
1629 :
1630 28101 : 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 20 : _fe_problem.getAuxiliarySystem().solution().close();
1636 20 : _fe_problem.getAuxiliarySystem().system().update();
1637 20 : solutionOld().close();
1638 : }
1639 :
1640 : // We may have additional tagged vectors that also need to be accumulated
1641 28101 : _fe_problem.addCachedResidual(0);
1642 28101 : }
1643 :
1644 : void
1645 4497 : NonlinearSystemBase::overwriteNodeFace(NumericVector<Number> & soln)
1646 : {
1647 : // Overwrite results from integrator in case we have explicit dynamics contact constraints
1648 4497 : auto & subproblem = _fe_problem.getDisplacedProblem()
1649 6898 : ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1650 6898 : : static_cast<SubProblem &>(_fe_problem);
1651 4497 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1652 :
1653 6898 : for (const auto & it : penetration_locators)
1654 : {
1655 2401 : PenetrationLocator & pen_loc = *(it.second);
1656 :
1657 2401 : const auto & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1658 2401 : const BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1659 2401 : const BoundaryID primary_boundary = pen_loc._primary_boundary;
1660 :
1661 2401 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, true))
1662 : {
1663 : const auto & constraints =
1664 2401 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, true);
1665 39217 : for (const auto i : index_range(secondary_nodes))
1666 : {
1667 36816 : const auto secondary_node_num = secondary_nodes[i];
1668 36816 : const Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1669 :
1670 36816 : if (secondary_node.processor_id() == processor_id())
1671 28816 : if (pen_loc._penetration_info[secondary_node_num])
1672 57632 : for (const auto & nfc : constraints)
1673 : {
1674 28816 : 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 57632 : if (nfc->secondaryBoundary() != secondary_boundary ||
1682 28816 : nfc->primaryBoundary() != primary_boundary)
1683 0 : continue;
1684 :
1685 28816 : nfc->overwriteBoundaryVariables(soln, secondary_node);
1686 : }
1687 : }
1688 : }
1689 : }
1690 4497 : soln.close();
1691 4497 : }
1692 :
1693 : void
1694 3374022 : NonlinearSystemBase::residualSetup()
1695 : {
1696 3374022 : TIME_SECTION("residualSetup", 3);
1697 :
1698 3374022 : SolverSystem::residualSetup();
1699 :
1700 7050842 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1701 : {
1702 3676820 : _kernels.residualSetup(tid);
1703 3676820 : _nodal_kernels.residualSetup(tid);
1704 3676820 : _dirac_kernels.residualSetup(tid);
1705 3676820 : if (_doing_dg)
1706 35297 : _dg_kernels.residualSetup(tid);
1707 3676820 : _interface_kernels.residualSetup(tid);
1708 3676820 : _element_dampers.residualSetup(tid);
1709 3676820 : _nodal_dampers.residualSetup(tid);
1710 3676820 : _integrated_bcs.residualSetup(tid);
1711 : }
1712 3374022 : _scalar_kernels.residualSetup();
1713 3374022 : _constraints.residualSetup();
1714 3374022 : _general_dampers.residualSetup();
1715 3374022 : _nodal_bcs.residualSetup();
1716 :
1717 : // Avoid recursion
1718 3374022 : if (this == &_fe_problem.currentNonlinearSystem())
1719 3281482 : _fe_problem.residualSetup();
1720 3374022 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
1721 3374022 : }
1722 :
1723 : void
1724 3278015 : NonlinearSystemBase::computeResidualInternal(const std::set<TagID> & tags)
1725 : {
1726 : parallel_object_only();
1727 :
1728 3278015 : TIME_SECTION("computeResidualInternal", 3);
1729 :
1730 3278015 : residualSetup();
1731 :
1732 3278015 : 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 3278015 : std::vector<UserObject *> uos;
1737 3278015 : _fe_problem.theWarehouse()
1738 6556030 : .query()
1739 3278015 : .condition<AttribSystem>("UserObject")
1740 3278015 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1741 3278015 : .queryInto(uos);
1742 3278015 : for (auto & uo : uos)
1743 0 : uo->residualSetup();
1744 3278015 : for (auto & uo : uos)
1745 : {
1746 0 : uo->initialize();
1747 0 : uo->execute();
1748 0 : uo->finalize();
1749 : }
1750 :
1751 : // reinit scalar variables
1752 6850609 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1753 3572594 : _fe_problem.reinitScalars(tid);
1754 :
1755 : // residual contributions from the domain
1756 : PARALLEL_TRY
1757 : {
1758 3278015 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1759 :
1760 3278015 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
1761 :
1762 3278015 : ComputeResidualThread cr(_fe_problem, tags);
1763 3278015 : 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 3277994 : if (_fe_problem.haveFV())
1771 : {
1772 : ComputeFVFluxResidualThread<FVRange> fvr(
1773 149562 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
1774 149562 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
1775 149562 : Threads::parallel_reduce(faces, fvr);
1776 149562 : }
1777 3277986 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
1778 3277986 : 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 3277986 : }
1786 :
1787 3277986 : unsigned int n_threads = libMesh::n_threads();
1788 6850546 : for (unsigned int i = 0; i < n_threads;
1789 : i++) // Add any cached residuals that might be hanging around
1790 3572560 : _fe_problem.addCachedResidual(i);
1791 3277986 : }
1792 3277986 : 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 3277865 : if (_scalar_kernels.hasActiveObjects())
1799 : {
1800 73607 : 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 73607 : if (!tags.size() || tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL))
1807 73586 : scalar_kernel_warehouse = &_scalar_kernels;
1808 21 : else if (tags.size() == 1)
1809 : scalar_kernel_warehouse =
1810 12 : &(_scalar_kernels.getVectorTagObjectWarehouse(*(tags.begin()), 0));
1811 : else
1812 : // scalar_kernels is not threading
1813 9 : scalar_kernel_warehouse = &(_scalar_kernels.getVectorTagsObjectWarehouse(tags, 0));
1814 :
1815 73607 : bool have_scalar_contributions = false;
1816 73607 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
1817 270724 : for (const auto & scalar_kernel : scalars)
1818 : {
1819 197117 : scalar_kernel->reinit();
1820 197117 : const std::vector<dof_id_type> & dof_indices = scalar_kernel->variable().dofIndices();
1821 197117 : const DofMap & dof_map = scalar_kernel->variable().dofMap();
1822 197117 : const dof_id_type first_dof = dof_map.first_dof();
1823 197117 : const dof_id_type end_dof = dof_map.end_dof();
1824 227550 : for (dof_id_type dof : dof_indices)
1825 : {
1826 197694 : if (dof >= first_dof && dof < end_dof)
1827 : {
1828 167261 : scalar_kernel->computeResidual();
1829 167261 : have_scalar_contributions = true;
1830 167261 : break;
1831 : }
1832 : }
1833 : }
1834 73607 : if (have_scalar_contributions)
1835 64520 : _fe_problem.addResidualScalar();
1836 73607 : }
1837 : }
1838 3277865 : PARALLEL_CATCH;
1839 :
1840 : // residual contributions from Block NodalKernels
1841 : PARALLEL_TRY
1842 : {
1843 3277865 : if (_nodal_kernels.hasActiveBlockObjects())
1844 : {
1845 15596 : TIME_SECTION("NodalKernels", 3 /*, "Computing NodalKernels"*/);
1846 :
1847 15596 : ComputeNodalKernelsThread cnk(_fe_problem, _nodal_kernels, tags);
1848 :
1849 15596 : const ConstNodeRange & range = _fe_problem.getCurrentAlgebraicNodeRange();
1850 :
1851 15596 : if (range.begin() != range.end())
1852 : {
1853 15596 : _fe_problem.reinitNode(*range.begin(), 0);
1854 :
1855 15596 : Threads::parallel_reduce(range, cnk);
1856 :
1857 15596 : unsigned int n_threads = libMesh::n_threads();
1858 38526 : for (unsigned int i = 0; i < n_threads;
1859 : i++) // Add any cached residuals that might be hanging around
1860 22930 : _fe_problem.addCachedResidual(i);
1861 : }
1862 15596 : }
1863 : }
1864 3277865 : PARALLEL_CATCH;
1865 :
1866 3277865 : 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 3277815 : if (_nodal_kernels.hasActiveBoundaryObjects())
1876 : {
1877 1004 : TIME_SECTION("NodalKernelBCs", 3 /*, "Computing NodalKernelBCs"*/);
1878 :
1879 1004 : ComputeNodalKernelBcsThread cnk(_fe_problem, _nodal_kernels, tags);
1880 :
1881 1004 : const ConstBndNodeRange & bnd_node_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
1882 :
1883 1004 : Threads::parallel_reduce(bnd_node_range, cnk);
1884 :
1885 1004 : unsigned int n_threads = libMesh::n_threads();
1886 2094 : for (unsigned int i = 0; i < n_threads;
1887 : i++) // Add any cached residuals that might be hanging around
1888 1090 : _fe_problem.addCachedResidual(i);
1889 1004 : }
1890 : }
1891 3277815 : PARALLEL_CATCH;
1892 :
1893 3277815 : mortarConstraints(Moose::ComputeType::Residual, tags, {});
1894 :
1895 3277815 : if (_residual_copy.get())
1896 : {
1897 0 : _Re_non_time->close();
1898 0 : _Re_non_time->localize(*_residual_copy);
1899 : }
1900 :
1901 3277815 : if (_need_residual_ghosted)
1902 : {
1903 13234 : _Re_non_time->close();
1904 13234 : *_residual_ghosted = *_Re_non_time;
1905 13234 : _residual_ghosted->close();
1906 : }
1907 :
1908 3277815 : PARALLEL_TRY { computeDiracContributions(tags, false); }
1909 3277803 : PARALLEL_CATCH;
1910 :
1911 3277803 : if (_fe_problem._has_constraints)
1912 : {
1913 23285 : PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1914 23285 : PARALLEL_CATCH;
1915 23285 : _Re_non_time->close();
1916 : }
1917 :
1918 : // Add in Residual contributions from other Constraints
1919 3277803 : if (_fe_problem._has_constraints)
1920 : {
1921 : PARALLEL_TRY
1922 : {
1923 : // Undisplaced Constraints
1924 23285 : constraintResiduals(*_Re_non_time, false);
1925 :
1926 : // Displaced Constraints
1927 23285 : if (_fe_problem.getDisplacedProblem())
1928 4816 : constraintResiduals(*_Re_non_time, true);
1929 :
1930 23285 : if (_fe_problem.computingNonlinearResid())
1931 10581 : _constraints.residualEnd();
1932 : }
1933 23285 : PARALLEL_CATCH;
1934 23285 : _Re_non_time->close();
1935 : }
1936 :
1937 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
1938 : // counters
1939 3277803 : _app.solutionInvalidity().syncIteration();
1940 3277803 : _app.solutionInvalidity().solutionInvalidAccumulation();
1941 3278316 : }
1942 :
1943 : void
1944 3467 : NonlinearSystemBase::computeResidualAndJacobianInternal(const std::set<TagID> & vector_tags,
1945 : const std::set<TagID> & matrix_tags)
1946 : {
1947 3467 : TIME_SECTION("computeResidualAndJacobianInternal", 3);
1948 :
1949 : // Make matrix ready to use
1950 3467 : activeAllMatrixTags();
1951 :
1952 10401 : for (auto tag : matrix_tags)
1953 : {
1954 6934 : if (!hasMatrix(tag))
1955 3467 : continue;
1956 :
1957 3467 : auto & jacobian = getMatrix(tag);
1958 : // Necessary for speed
1959 3467 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
1960 : {
1961 3467 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
1962 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
1963 : PETSC_TRUE));
1964 3467 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
1965 0 : LibmeshPetscCall(
1966 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1967 3467 : 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 3467 : 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 3467 : std::vector<UserObject *> uos;
1979 3467 : _fe_problem.theWarehouse()
1980 6934 : .query()
1981 3467 : .condition<AttribSystem>("UserObject")
1982 3467 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1983 3467 : .queryInto(uos);
1984 3467 : for (auto & uo : uos)
1985 0 : uo->residualSetup();
1986 3467 : for (auto & uo : uos)
1987 : {
1988 0 : uo->initialize();
1989 0 : uo->execute();
1990 0 : uo->finalize();
1991 : }
1992 :
1993 : // reinit scalar variables
1994 7273 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1995 3806 : _fe_problem.reinitScalars(tid);
1996 :
1997 : // residual contributions from the domain
1998 : PARALLEL_TRY
1999 : {
2000 3467 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
2001 :
2002 3467 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
2003 :
2004 3467 : ComputeResidualAndJacobianThread crj(_fe_problem, vector_tags, matrix_tags);
2005 3467 : Threads::parallel_reduce(elem_range, crj);
2006 :
2007 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
2008 3467 : if (_fe_problem.haveFV())
2009 : {
2010 : ComputeFVFluxRJThread<FVRange> fvrj(
2011 1451 : _fe_problem, this->number(), vector_tags, matrix_tags, /*on_displaced=*/false);
2012 1451 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
2013 1451 : Threads::parallel_reduce(faces, fvrj);
2014 1451 : }
2015 3467 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
2016 3467 : 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 3467 : }
2024 :
2025 3467 : mortarConstraints(Moose::ComputeType::ResidualAndJacobian, vector_tags, matrix_tags);
2026 :
2027 3467 : unsigned int n_threads = libMesh::n_threads();
2028 7273 : for (unsigned int i = 0; i < n_threads;
2029 : i++) // Add any cached residuals that might be hanging around
2030 : {
2031 3806 : _fe_problem.addCachedResidual(i);
2032 3806 : _fe_problem.addCachedJacobian(i);
2033 : }
2034 3467 : }
2035 3467 : PARALLEL_CATCH;
2036 3467 : }
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 3277803 : 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 3277803 : if (_has_save_in)
2068 309 : _fe_problem.getAuxiliarySystem().solution().close();
2069 :
2070 : // Select nodal kernels
2071 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2072 :
2073 3277803 : if (tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL) || !tags.size())
2074 3239752 : nbc_warehouse = &_nodal_bcs;
2075 38051 : else if (tags.size() == 1)
2076 18286 : nbc_warehouse = &(_nodal_bcs.getVectorTagObjectWarehouse(*(tags.begin()), 0));
2077 : else
2078 19765 : nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0));
2079 :
2080 : // Return early if there is no nodal kernel
2081 3277803 : if (!nbc_warehouse->size())
2082 418442 : return;
2083 :
2084 : PARALLEL_TRY
2085 : {
2086 2859361 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2087 :
2088 2859361 : if (!bnd_nodes.empty())
2089 : {
2090 2859050 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2091 :
2092 150005303 : for (const auto & bnode : bnd_nodes)
2093 : {
2094 147146253 : BoundaryID boundary_id = bnode->_bnd_id;
2095 147146253 : Node * node = bnode->_node;
2096 :
2097 260318635 : if (node->processor_id() == processor_id() &&
2098 113172382 : nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
2099 : {
2100 : // reinit variables in nodes
2101 57916007 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2102 :
2103 57916007 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2104 121298770 : for (const auto & nbc : bcs)
2105 63382763 : if (nbc->shouldApply())
2106 63381359 : nbc->computeResidual();
2107 : }
2108 : }
2109 2859050 : }
2110 : }
2111 2859361 : PARALLEL_CATCH;
2112 :
2113 2859361 : if (_Re_time)
2114 2536777 : _Re_time->close();
2115 2859361 : _Re_non_time->close();
2116 : }
2117 :
2118 : void
2119 3467 : NonlinearSystemBase::computeNodalBCsResidualAndJacobian()
2120 : {
2121 : PARALLEL_TRY
2122 : {
2123 3467 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2124 :
2125 3467 : if (!bnd_nodes.empty())
2126 : {
2127 3467 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2128 :
2129 86555 : for (const auto & bnode : bnd_nodes)
2130 : {
2131 83088 : BoundaryID boundary_id = bnode->_bnd_id;
2132 83088 : Node * node = bnode->_node;
2133 :
2134 83088 : if (node->processor_id() == processor_id())
2135 : {
2136 : // reinit variables in nodes
2137 62704 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2138 62704 : if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
2139 : {
2140 12501 : const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
2141 25002 : for (const auto & nbc : bcs)
2142 12501 : if (nbc->shouldApply())
2143 12501 : nbc->computeResidualAndJacobian();
2144 : }
2145 : }
2146 : }
2147 3467 : }
2148 : }
2149 3467 : PARALLEL_CATCH;
2150 :
2151 : // Set the cached NodalBCBase values in the Jacobian matrix
2152 3467 : _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{});
2153 3467 : }
2154 :
2155 : void
2156 717 : NonlinearSystemBase::getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs)
2157 : {
2158 717 : const Node & node = _mesh.nodeRef(node_id);
2159 717 : unsigned int s = number();
2160 717 : if (node.has_dofs(s))
2161 : {
2162 1472 : for (unsigned int v = 0; v < nVariables(); v++)
2163 1472 : for (unsigned int c = 0; c < node.n_comp(s, v); c++)
2164 717 : dofs.push_back(node.dof_number(s, v, c));
2165 : }
2166 717 : }
2167 :
2168 : void
2169 5033 : NonlinearSystemBase::findImplicitGeometricCouplingEntries(
2170 : GeometricSearchData & geom_search_data,
2171 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & graph)
2172 : {
2173 5033 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
2174 5033 : const auto & nearest_node_locators = geom_search_data._nearest_node_locators;
2175 7162 : for (const auto & it : nearest_node_locators)
2176 : {
2177 2129 : std::vector<dof_id_type> & secondary_nodes = it.second->_secondary_nodes;
2178 :
2179 10915 : for (const auto & secondary_node : secondary_nodes)
2180 : {
2181 8786 : std::set<dof_id_type> unique_secondary_indices;
2182 8786 : std::set<dof_id_type> unique_primary_indices;
2183 :
2184 8786 : auto node_to_elem_pair = node_to_elem_map.find(secondary_node);
2185 8786 : if (node_to_elem_pair != node_to_elem_map.end())
2186 : {
2187 8666 : 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 23314 : for (const auto & cur_elem : elems)
2191 : {
2192 14648 : std::vector<dof_id_type> dof_indices;
2193 14648 : dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2194 :
2195 73240 : for (const auto & dof : dof_indices)
2196 58592 : unique_secondary_indices.insert(dof);
2197 14648 : }
2198 : }
2199 :
2200 8786 : std::vector<dof_id_type> primary_nodes = it.second->_neighbor_nodes[secondary_node];
2201 :
2202 49576 : for (const auto & primary_node : primary_nodes)
2203 : {
2204 40790 : 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 40790 : 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 107078 : for (const auto & cur_elem : primary_node_elems)
2211 : {
2212 66288 : std::vector<dof_id_type> dof_indices;
2213 66288 : dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2214 :
2215 331440 : for (const auto & dof : dof_indices)
2216 265152 : unique_primary_indices.insert(dof);
2217 66288 : }
2218 : }
2219 :
2220 55414 : for (const auto & secondary_id : unique_secondary_indices)
2221 506684 : for (const auto & primary_id : unique_primary_indices)
2222 : {
2223 460056 : graph[secondary_id].push_back(primary_id);
2224 460056 : graph[primary_id].push_back(secondary_id);
2225 : }
2226 8786 : }
2227 : }
2228 :
2229 : // handle node-to-node constraints
2230 5033 : const auto & ncs = _constraints.getActiveNodalConstraints();
2231 5190 : for (const auto & nc : ncs)
2232 : {
2233 157 : std::vector<dof_id_type> primary_dofs;
2234 157 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
2235 314 : for (const auto & node_id : primary_node_ids)
2236 : {
2237 157 : Node * node = _mesh.queryNodePtr(node_id);
2238 157 : if (node && node->processor_id() == this->processor_id())
2239 : {
2240 142 : getNodeDofs(node_id, primary_dofs);
2241 : }
2242 : }
2243 :
2244 157 : _communicator.allgather(primary_dofs);
2245 :
2246 157 : std::vector<dof_id_type> secondary_dofs;
2247 157 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
2248 732 : for (const auto & node_id : secondary_node_ids)
2249 : {
2250 575 : Node * node = _mesh.queryNodePtr(node_id);
2251 575 : if (node && node->processor_id() == this->processor_id())
2252 : {
2253 575 : getNodeDofs(node_id, secondary_dofs);
2254 : }
2255 : }
2256 :
2257 157 : _communicator.allgather(secondary_dofs);
2258 :
2259 314 : for (const auto & primary_id : primary_dofs)
2260 870 : for (const auto & secondary_id : secondary_dofs)
2261 : {
2262 713 : graph[primary_id].push_back(secondary_id);
2263 713 : graph[secondary_id].push_back(primary_id);
2264 : }
2265 157 : }
2266 :
2267 : // Make every entry sorted and unique
2268 29785 : for (auto & it : graph)
2269 : {
2270 24752 : std::vector<dof_id_type> & row = it.second;
2271 24752 : std::sort(row.begin(), row.end());
2272 24752 : std::vector<dof_id_type>::iterator uit = std::unique(row.begin(), row.end());
2273 24752 : row.resize(uit - row.begin());
2274 : }
2275 5033 : }
2276 :
2277 : void
2278 4096 : NonlinearSystemBase::addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data)
2279 : {
2280 4096 : 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 4096 : auto & jacobian = getMatrix(systemMatrixTag());
2286 :
2287 4096 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
2288 :
2289 4096 : findImplicitGeometricCouplingEntries(geom_search_data, graph);
2290 :
2291 27971 : for (const auto & it : graph)
2292 : {
2293 23875 : dof_id_type dof = it.first;
2294 23875 : const auto & row = it.second;
2295 :
2296 281261 : for (const auto & coupled_dof : row)
2297 257386 : jacobian.add(dof, coupled_dof, 0);
2298 : }
2299 4096 : }
2300 :
2301 : void
2302 4688 : NonlinearSystemBase::constraintJacobians(const SparseMatrix<Number> & jacobian_to_view,
2303 : bool displaced)
2304 : {
2305 4688 : if (!hasMatrix(systemMatrixTag()))
2306 0 : mooseError("A system matrix is required");
2307 :
2308 4688 : auto & jacobian = getMatrix(systemMatrixTag());
2309 :
2310 4688 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2311 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2312 : MAT_NEW_NONZERO_ALLOCATION_ERR,
2313 : PETSC_FALSE));
2314 4688 : if (_fe_problem.ignoreZerosInJacobian())
2315 0 : LibmeshPetscCall(MatSetOption(
2316 : static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2317 :
2318 4688 : 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 2142 : auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
2325 5759 : : static_cast<SubProblem &>(_fe_problem);
2326 4688 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
2327 :
2328 : bool constraints_applied;
2329 4688 : if (!_assemble_constraints_separately)
2330 4688 : constraints_applied = false;
2331 6754 : for (const auto & it : penetration_locators)
2332 : {
2333 2066 : 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 2066 : PenetrationLocator & pen_loc = *(it.second);
2340 :
2341 2066 : std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
2342 :
2343 2066 : BoundaryID secondary_boundary = pen_loc._secondary_boundary;
2344 2066 : BoundaryID primary_boundary = pen_loc._primary_boundary;
2345 :
2346 2066 : zero_rows.clear();
2347 2066 : if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
2348 : {
2349 : const auto & constraints =
2350 1128 : _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
2351 :
2352 5726 : for (const auto & secondary_node_num : secondary_nodes)
2353 : {
2354 4598 : Node & secondary_node = _mesh.nodeRef(secondary_node_num);
2355 :
2356 4598 : if (secondary_node.processor_id() == processor_id())
2357 : {
2358 4184 : if (pen_loc._penetration_info[secondary_node_num])
2359 : {
2360 4184 : PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
2361 :
2362 4184 : reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
2363 4184 : _fe_problem.reinitOffDiagScalars(0);
2364 :
2365 8368 : for (const auto & nfc : constraints)
2366 : {
2367 4184 : 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 8368 : if (nfc->secondaryBoundary() != secondary_boundary ||
2374 4184 : nfc->primaryBoundary() != primary_boundary)
2375 0 : continue;
2376 :
2377 4184 : nfc->_jacobian = &jacobian_to_view;
2378 :
2379 4184 : if (nfc->shouldApply())
2380 : {
2381 4184 : constraints_applied = true;
2382 :
2383 4184 : nfc->prepareShapes(nfc->variable().number());
2384 4184 : nfc->prepareNeighborShapes(nfc->variable().number());
2385 :
2386 4184 : nfc->computeJacobian();
2387 :
2388 4184 : if (nfc->overwriteSecondaryJacobian())
2389 : {
2390 : // Add this variable's dof's row to be zeroed
2391 4184 : zero_rows.push_back(nfc->variable().nodalDofIndex());
2392 : }
2393 :
2394 4184 : 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 4184 : nfc->overwriteSecondaryJacobian() ? 1. : nfc->variable().scalingFactor();
2402 :
2403 : // Cache the jacobian block for the secondary side
2404 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2405 4184 : nfc->_Kee,
2406 : secondary_dofs,
2407 4184 : nfc->_connected_dof_indices,
2408 : scaling_factor);
2409 :
2410 : // Cache Ken, Kne, Knn
2411 4184 : 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 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2416 4184 : nfc->_Ken,
2417 : secondary_dofs,
2418 4184 : nfc->primaryVariable().dofIndicesNeighbor(),
2419 : scaling_factor);
2420 :
2421 : // Use _connected_dof_indices to get all the correct columns
2422 8368 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2423 4184 : nfc->_Kne,
2424 4184 : nfc->primaryVariable().dofIndicesNeighbor(),
2425 4184 : nfc->_connected_dof_indices,
2426 4184 : nfc->primaryVariable().scalingFactor());
2427 :
2428 : // We've handled Ken and Kne, finally handle Knn
2429 4184 : _fe_problem.cacheJacobianNeighbor(0);
2430 : }
2431 :
2432 : // Do the off-diagonals next
2433 4184 : const std::vector<MooseVariableFEBase *> coupled_vars = nfc->getCoupledMooseVars();
2434 8368 : for (const auto & jvar : coupled_vars)
2435 : {
2436 : // Only compute jacobians for nonlinear variables
2437 4184 : 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 4256 : if (nfc->variable().number() == jvar->number() ||
2443 144 : !_fe_problem.areCoupled(
2444 72 : nfc->variable().number(), jvar->number(), this->number()))
2445 4112 : continue;
2446 :
2447 : // Need to zero out the matrices first
2448 72 : _fe_problem.prepareAssembly(0);
2449 :
2450 72 : nfc->prepareShapes(nfc->variable().number());
2451 72 : nfc->prepareNeighborShapes(jvar->number());
2452 :
2453 72 : nfc->computeOffDiagJacobian(jvar->number());
2454 :
2455 : // Cache the jacobian block for the secondary side
2456 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2457 72 : nfc->_Kee,
2458 : secondary_dofs,
2459 72 : nfc->_connected_dof_indices,
2460 : scaling_factor);
2461 :
2462 : // Cache Ken, Kne, Knn
2463 72 : 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 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2468 72 : nfc->_Ken,
2469 : secondary_dofs,
2470 72 : jvar->dofIndicesNeighbor(),
2471 : scaling_factor);
2472 :
2473 : // Use _connected_dof_indices to get all the correct columns
2474 144 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2475 72 : nfc->_Kne,
2476 72 : nfc->variable().dofIndicesNeighbor(),
2477 72 : nfc->_connected_dof_indices,
2478 72 : nfc->variable().scalingFactor());
2479 :
2480 : // We've handled Ken and Kne, finally handle Knn
2481 72 : _fe_problem.cacheJacobianNeighbor(0);
2482 : }
2483 : }
2484 4184 : }
2485 : }
2486 : }
2487 : }
2488 : }
2489 : }
2490 2066 : 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 4688 : if (!_assemble_constraints_separately)
2510 : {
2511 : // See if constraints were applied anywhere
2512 4688 : _communicator.max(constraints_applied);
2513 :
2514 4688 : if (constraints_applied)
2515 : {
2516 1054 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2517 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2518 : PETSC_TRUE));
2519 :
2520 1054 : jacobian.close();
2521 1054 : jacobian.zero_rows(zero_rows, 0.0);
2522 1054 : jacobian.close();
2523 1054 : _fe_problem.addCachedJacobian(0);
2524 1054 : jacobian.close();
2525 : }
2526 : }
2527 :
2528 4688 : THREAD_ID tid = 0;
2529 : // go over element-element constraint interface
2530 4688 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
2531 4688 : 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 4688 : std::set<dof_id_type> unique_secondary_node_ids;
2577 4688 : constraints_applied = false;
2578 19936 : for (const auto & secondary_id : _mesh.meshSubdomains())
2579 : {
2580 74132 : for (const auto & primary_id : _mesh.meshSubdomains())
2581 : {
2582 58884 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
2583 : {
2584 : const auto & constraints =
2585 180 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
2586 :
2587 : // get unique set of ids of all nodes on current block
2588 180 : unique_secondary_node_ids.clear();
2589 180 : const MeshBase & meshhelper = _mesh.getMesh();
2590 360 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
2591 17340 : meshhelper.active_subdomain_elements_end(secondary_id)))
2592 : {
2593 49140 : for (auto & n : elem->node_ref_range())
2594 40740 : unique_secondary_node_ids.insert(n.id());
2595 180 : }
2596 :
2597 14040 : for (auto secondary_node_id : unique_secondary_node_ids)
2598 : {
2599 13860 : const Node & secondary_node = _mesh.nodeRef(secondary_node_id);
2600 : // check if secondary node is on current processor
2601 13860 : if (secondary_node.processor_id() == processor_id())
2602 : {
2603 : // This reinits the variables that exist on the secondary node
2604 11088 : _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 11088 : _fe_problem.prepareAssembly(0);
2609 11088 : _fe_problem.reinitOffDiagScalars(0);
2610 :
2611 23584 : for (const auto & nec : constraints)
2612 : {
2613 12496 : if (nec->shouldApply())
2614 : {
2615 5680 : constraints_applied = true;
2616 :
2617 5680 : nec->_jacobian = &jacobian_to_view;
2618 5680 : nec->prepareShapes(nec->variable().number());
2619 5680 : nec->prepareNeighborShapes(nec->variable().number());
2620 :
2621 5680 : nec->computeJacobian();
2622 :
2623 5680 : 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 5680 : std::vector<dof_id_type> secondary_dofs(1, nec->variable().nodalDofIndex());
2630 :
2631 : // Cache the jacobian block for the secondary side
2632 11360 : nec->addJacobian(_fe_problem.assembly(0, number()),
2633 5680 : nec->_Kee,
2634 : secondary_dofs,
2635 5680 : nec->_connected_dof_indices,
2636 5680 : nec->variable().scalingFactor());
2637 :
2638 : // Cache the jacobian block for the primary side
2639 11360 : nec->addJacobian(_fe_problem.assembly(0, number()),
2640 5680 : nec->_Kne,
2641 5680 : nec->primaryVariable().dofIndicesNeighbor(),
2642 5680 : nec->_connected_dof_indices,
2643 5680 : nec->primaryVariable().scalingFactor());
2644 :
2645 5680 : _fe_problem.cacheJacobian(0);
2646 5680 : _fe_problem.cacheJacobianNeighbor(0);
2647 :
2648 : // Do the off-diagonals next
2649 5680 : const std::vector<MooseVariableFEBase *> coupled_vars = nec->getCoupledMooseVars();
2650 11440 : for (const auto & jvar : coupled_vars)
2651 : {
2652 : // Only compute jacobians for nonlinear variables
2653 5760 : if (jvar->kind() != Moose::VAR_SOLVER)
2654 80 : continue;
2655 :
2656 : // Only compute Jacobian entries if this coupling is being used by the
2657 : // preconditioner
2658 6320 : if (nec->variable().number() == jvar->number() ||
2659 1280 : !_fe_problem.areCoupled(
2660 640 : nec->variable().number(), jvar->number(), this->number()))
2661 5040 : continue;
2662 :
2663 : // Need to zero out the matrices first
2664 640 : _fe_problem.prepareAssembly(0);
2665 :
2666 640 : nec->prepareShapes(nec->variable().number());
2667 640 : nec->prepareNeighborShapes(jvar->number());
2668 :
2669 640 : nec->computeOffDiagJacobian(jvar->number());
2670 :
2671 : // Cache the jacobian block for the secondary side
2672 1280 : nec->addJacobian(_fe_problem.assembly(0, number()),
2673 640 : nec->_Kee,
2674 : secondary_dofs,
2675 640 : nec->_connected_dof_indices,
2676 640 : nec->variable().scalingFactor());
2677 :
2678 : // Cache the jacobian block for the primary side
2679 1280 : nec->addJacobian(_fe_problem.assembly(0, number()),
2680 640 : nec->_Kne,
2681 640 : nec->variable().dofIndicesNeighbor(),
2682 640 : nec->_connected_dof_indices,
2683 640 : nec->variable().scalingFactor());
2684 :
2685 640 : _fe_problem.cacheJacobian(0);
2686 640 : _fe_problem.cacheJacobianNeighbor(0);
2687 : }
2688 5680 : }
2689 : }
2690 : }
2691 : }
2692 : }
2693 : }
2694 : }
2695 : // See if constraints were applied anywhere
2696 4688 : _communicator.max(constraints_applied);
2697 :
2698 4688 : if (constraints_applied)
2699 : {
2700 180 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2701 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2702 : PETSC_TRUE));
2703 :
2704 180 : jacobian.close();
2705 180 : jacobian.zero_rows(zero_rows, 0.0);
2706 180 : jacobian.close();
2707 180 : _fe_problem.addCachedJacobian(0);
2708 180 : jacobian.close();
2709 : }
2710 4688 : }
2711 :
2712 : void
2713 517237 : NonlinearSystemBase::computeScalarKernelsJacobians(const std::set<TagID> & tags)
2714 : {
2715 : MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
2716 :
2717 517237 : if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2718 507666 : scalar_kernel_warehouse = &_scalar_kernels;
2719 9571 : else if (tags.size() == 1)
2720 8319 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2721 : else
2722 1252 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagsObjectWarehouse(tags, 0));
2723 :
2724 : // Compute the diagonal block for scalar variables
2725 517237 : if (scalar_kernel_warehouse->hasActiveObjects())
2726 : {
2727 17517 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
2728 :
2729 17517 : _fe_problem.reinitScalars(/*tid=*/0);
2730 :
2731 17517 : _fe_problem.reinitOffDiagScalars(/*_tid*/ 0);
2732 :
2733 17517 : bool have_scalar_contributions = false;
2734 62793 : for (const auto & kernel : scalars)
2735 : {
2736 45276 : kernel->reinit();
2737 45276 : const std::vector<dof_id_type> & dof_indices = kernel->variable().dofIndices();
2738 45276 : const DofMap & dof_map = kernel->variable().dofMap();
2739 45276 : const dof_id_type first_dof = dof_map.first_dof();
2740 45276 : const dof_id_type end_dof = dof_map.end_dof();
2741 51417 : for (dof_id_type dof : dof_indices)
2742 : {
2743 45318 : if (dof >= first_dof && dof < end_dof)
2744 : {
2745 39177 : kernel->computeJacobian();
2746 39177 : _fe_problem.addJacobianOffDiagScalar(kernel->variable().number());
2747 39177 : have_scalar_contributions = true;
2748 39177 : break;
2749 : }
2750 : }
2751 : }
2752 :
2753 17517 : if (have_scalar_contributions)
2754 15594 : _fe_problem.addJacobianScalar();
2755 : }
2756 517237 : }
2757 :
2758 : void
2759 533297 : NonlinearSystemBase::jacobianSetup()
2760 : {
2761 533297 : SolverSystem::jacobianSetup();
2762 :
2763 1115850 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2764 : {
2765 582553 : _kernels.jacobianSetup(tid);
2766 582553 : _nodal_kernels.jacobianSetup(tid);
2767 582553 : _dirac_kernels.jacobianSetup(tid);
2768 582553 : if (_doing_dg)
2769 2948 : _dg_kernels.jacobianSetup(tid);
2770 582553 : _interface_kernels.jacobianSetup(tid);
2771 582553 : _element_dampers.jacobianSetup(tid);
2772 582553 : _nodal_dampers.jacobianSetup(tid);
2773 582553 : _integrated_bcs.jacobianSetup(tid);
2774 : }
2775 533297 : _scalar_kernels.jacobianSetup();
2776 533297 : _constraints.jacobianSetup();
2777 533297 : _general_dampers.jacobianSetup();
2778 533297 : _nodal_bcs.jacobianSetup();
2779 :
2780 : // Avoid recursion
2781 533297 : if (this == &_fe_problem.currentNonlinearSystem())
2782 517237 : _fe_problem.jacobianSetup();
2783 533297 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
2784 533297 : }
2785 :
2786 : void
2787 517237 : NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
2788 : {
2789 517237 : TIME_SECTION("computeJacobianInternal", 3);
2790 :
2791 517237 : _fe_problem.setCurrentNonlinearSystem(number());
2792 :
2793 : // Make matrix ready to use
2794 517237 : activeAllMatrixTags();
2795 :
2796 1543884 : for (auto tag : tags)
2797 : {
2798 1026647 : if (!hasMatrix(tag))
2799 507666 : continue;
2800 :
2801 518981 : auto & jacobian = getMatrix(tag);
2802 : // Necessary for speed
2803 518981 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2804 : {
2805 517755 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2806 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2807 : PETSC_TRUE));
2808 517755 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2809 13371 : LibmeshPetscCall(
2810 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2811 517755 : 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 517237 : 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 517237 : std::vector<UserObject *> uos;
2823 517237 : _fe_problem.theWarehouse()
2824 1034474 : .query()
2825 517237 : .condition<AttribSystem>("UserObject")
2826 517237 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
2827 517237 : .queryInto(uos);
2828 517237 : for (auto & uo : uos)
2829 0 : uo->jacobianSetup();
2830 517237 : for (auto & uo : uos)
2831 : {
2832 0 : uo->initialize();
2833 0 : uo->execute();
2834 0 : uo->finalize();
2835 : }
2836 :
2837 : // reinit scalar variables
2838 1082357 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2839 565120 : _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 517237 : computeScalarKernelsJacobians(tags);
2847 :
2848 : // Block restricted Nodal Kernels
2849 517237 : if (_nodal_kernels.hasActiveBlockObjects())
2850 : {
2851 3998 : ComputeNodalKernelJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
2852 3998 : const ConstNodeRange & range = _fe_problem.getCurrentAlgebraicNodeRange();
2853 3998 : Threads::parallel_reduce(range, cnkjt);
2854 :
2855 3998 : unsigned int n_threads = libMesh::n_threads();
2856 9262 : for (unsigned int i = 0; i < n_threads;
2857 : i++) // Add any cached jacobians that might be hanging around
2858 5264 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
2859 3998 : }
2860 :
2861 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
2862 517237 : 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 37322 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
2868 37322 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
2869 37322 : Threads::parallel_reduce(faces, fvj);
2870 37322 : }
2871 517237 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
2872 517237 : 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 517237 : }
2880 :
2881 517237 : mortarConstraints(Moose::ComputeType::Jacobian, {}, tags);
2882 :
2883 : // Get our element range for looping over
2884 517237 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
2885 :
2886 517237 : 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 602 : ComputeJacobianForScalingThread cj(_fe_problem, tags);
2894 602 : Threads::parallel_reduce(elem_range, cj);
2895 602 : unsigned int n_threads = libMesh::n_threads();
2896 1263 : for (unsigned int i = 0; i < n_threads;
2897 : i++) // Add any Jacobian contributions still hanging around
2898 661 : _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 602 : _fe_problem.checkExceptionAndStopSolve();
2905 :
2906 602 : closeTaggedMatrices(tags);
2907 :
2908 602 : return;
2909 602 : }
2910 :
2911 516635 : switch (_fe_problem.coupling())
2912 : {
2913 428711 : case Moose::COUPLING_DIAG:
2914 : {
2915 428711 : ComputeJacobianThread cj(_fe_problem, tags);
2916 428711 : Threads::parallel_reduce(elem_range, cj);
2917 :
2918 428703 : unsigned int n_threads = libMesh::n_threads();
2919 898990 : for (unsigned int i = 0; i < n_threads;
2920 : i++) // Add any Jacobian contributions still hanging around
2921 470287 : _fe_problem.addCachedJacobian(i);
2922 :
2923 : // Boundary restricted Nodal Kernels
2924 428703 : if (_nodal_kernels.hasActiveBoundaryObjects())
2925 : {
2926 46 : ComputeNodalKernelBCJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
2927 46 : const ConstBndNodeRange & bnd_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
2928 :
2929 46 : Threads::parallel_reduce(bnd_range, cnkjt);
2930 46 : unsigned int n_threads = libMesh::n_threads();
2931 96 : for (unsigned int i = 0; i < n_threads;
2932 : i++) // Add any cached jacobians that might be hanging around
2933 50 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
2934 46 : }
2935 428703 : }
2936 428703 : break;
2937 :
2938 87924 : default:
2939 : case Moose::COUPLING_CUSTOM:
2940 : {
2941 87924 : ComputeFullJacobianThread cj(_fe_problem, tags);
2942 87924 : Threads::parallel_reduce(elem_range, cj);
2943 87924 : unsigned int n_threads = libMesh::n_threads();
2944 :
2945 182083 : for (unsigned int i = 0; i < n_threads; i++)
2946 94163 : _fe_problem.addCachedJacobian(i);
2947 :
2948 : // Boundary restricted Nodal Kernels
2949 87920 : 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 87924 : }
2961 87920 : break;
2962 : }
2963 :
2964 516623 : 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 1032646 : if ((_fe_problem.restoreOriginalNonzeroPattern() || first) &&
2970 516027 : _add_implicit_geometric_coupling_entries_to_jacobian)
2971 : {
2972 3025 : first = false;
2973 3025 : addImplicitGeometricCouplingEntries(_fe_problem.geomSearchData());
2974 :
2975 3025 : if (_fe_problem.getDisplacedProblem())
2976 1071 : addImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData());
2977 : }
2978 : }
2979 516619 : 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 516517 : 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 3617 : auto & system_matrix = getMatrix(systemMatrixTag());
2991 : #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
2992 : SparseMatrix<Number> * view_jac_ptr;
2993 3617 : std::unique_ptr<SparseMatrix<Number>> hash_copy;
2994 3617 : if (system_matrix.use_hash_table())
2995 : {
2996 2309 : hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
2997 2309 : view_jac_ptr = hash_copy.get();
2998 : }
2999 : else
3000 1308 : view_jac_ptr = &system_matrix;
3001 3617 : auto & jacobian_to_view = *view_jac_ptr;
3002 : #else
3003 : auto & jacobian_to_view = system_matrix;
3004 : #endif
3005 3617 : if (&jacobian_to_view == &system_matrix)
3006 1308 : system_matrix.close();
3007 :
3008 : // Nodal Constraints
3009 3617 : enforceNodalConstraintsJacobian();
3010 :
3011 : // Undisplaced Constraints
3012 3617 : constraintJacobians(jacobian_to_view, false);
3013 :
3014 : // Displaced Constraints
3015 3617 : if (_fe_problem.getDisplacedProblem())
3016 1071 : constraintJacobians(jacobian_to_view, true);
3017 3617 : }
3018 : }
3019 516517 : 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 516517 : if (_has_diag_save_in)
3024 185 : _fe_problem.getAuxiliarySystem().solution().close();
3025 :
3026 : PARALLEL_TRY
3027 : {
3028 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
3029 : // Select nodal kernels
3030 516517 : if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
3031 506956 : nbc_warehouse = &_nodal_bcs;
3032 9561 : else if (tags.size() == 1)
3033 8309 : nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
3034 : else
3035 1252 : nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0));
3036 :
3037 516517 : 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 425223 : 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 425223 : std::map<std::string, std::set<unsigned int>> bc_involved_vars;
3050 425223 : const std::set<BoundaryID> & all_boundary_ids = _mesh.getBoundaryIDs();
3051 2136783 : 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 1711560 : if (nbc_warehouse->hasActiveBoundaryObjects(bid))
3056 : {
3057 841188 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
3058 1764069 : for (const auto & bc : bcs)
3059 : {
3060 : const std::vector<MooseVariableFEBase *> & coupled_moose_vars =
3061 922881 : 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 922881 : std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3066 980953 : for (const auto & coupled_var : coupled_moose_vars)
3067 58072 : if (coupled_var->kind() == Moose::VAR_SOLVER)
3068 56694 : var_set.insert(coupled_var->number());
3069 :
3070 922881 : 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 888607 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3079 463384 : _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 425223 : auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number());
3087 :
3088 : // Compute Jacobians for NodalBCBases
3089 425223 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3090 22285122 : for (const auto & bnode : bnd_nodes)
3091 : {
3092 21859899 : BoundaryID boundary_id = bnode->_bnd_id;
3093 21859899 : Node * node = bnode->_node;
3094 :
3095 32068359 : if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
3096 10208460 : node->processor_id() == processor_id())
3097 : {
3098 7961097 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3099 :
3100 7961097 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
3101 16919719 : for (const auto & bc : bcs)
3102 : {
3103 : // Get the set of involved MOOSE vars for this BC
3104 8958622 : 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 22482860 : for (const auto & it : coupling_entries)
3111 : {
3112 13524238 : 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 13524238 : if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
3119 9048006 : bc->computeOffDiagJacobian(jvar);
3120 : }
3121 :
3122 8958622 : const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars();
3123 8958993 : for (const auto & jvariable : coupled_scalar_vars)
3124 371 : if (hasScalarVariable(jvariable->name()))
3125 371 : bc->computeOffDiagJacobianScalar(jvariable->number());
3126 : }
3127 : }
3128 : } // end loop over boundary nodes
3129 :
3130 : // Set the cached NodalBCBase values in the Jacobian matrix
3131 425223 : _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{});
3132 425223 : }
3133 : }
3134 516517 : PARALLEL_CATCH;
3135 :
3136 516517 : 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 516517 : if (_has_nodalbc_diag_save_in)
3141 23 : _fe_problem.getAuxiliarySystem().solution().close();
3142 :
3143 516517 : if (hasDiagSaveIn())
3144 185 : _fe_problem.getAuxiliarySystem().update();
3145 :
3146 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
3147 : // counters
3148 516517 : _app.solutionInvalidity().syncIteration();
3149 516517 : _app.solutionInvalidity().solutionInvalidAccumulation();
3150 517933 : }
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 517237 : NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
3177 : {
3178 517237 : TIME_SECTION("computeJacobianTags", 5);
3179 :
3180 517237 : FloatingPointExceptionGuard fpe_guard(_app);
3181 :
3182 : try
3183 : {
3184 517237 : computeJacobianInternal(tags);
3185 : }
3186 106 : 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 102 : }
3192 517229 : }
3193 :
3194 : void
3195 228 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
3196 : {
3197 228 : _nl_matrix_tags.clear();
3198 :
3199 228 : auto & tags = _fe_problem.getMatrixTags();
3200 684 : for (auto & tag : tags)
3201 456 : _nl_matrix_tags.insert(tag.second);
3202 :
3203 228 : computeJacobianBlocks(blocks, _nl_matrix_tags);
3204 228 : }
3205 :
3206 : void
3207 588 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
3208 : const std::set<TagID> & tags)
3209 : {
3210 588 : TIME_SECTION("computeJacobianBlocks", 3);
3211 588 : FloatingPointExceptionGuard fpe_guard(_app);
3212 :
3213 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3214 : {
3215 1512 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3216 :
3217 1512 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3218 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
3219 : PETSC_TRUE));
3220 1512 : 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 1512 : jacobian.zero();
3226 : }
3227 :
3228 1233 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3229 645 : _fe_problem.reinitScalars(tid);
3230 :
3231 : PARALLEL_TRY
3232 : {
3233 588 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
3234 588 : ComputeJacobianBlocksThread cjb(_fe_problem, blocks, tags);
3235 588 : Threads::parallel_reduce(elem_range, cjb);
3236 588 : }
3237 588 : PARALLEL_CATCH;
3238 :
3239 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3240 1512 : blocks[i]->_jacobian.close();
3241 :
3242 2100 : for (unsigned int i = 0; i < blocks.size(); i++)
3243 : {
3244 1512 : libMesh::System & precond_system = blocks[i]->_precond_system;
3245 1512 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3246 :
3247 1512 : unsigned int ivar = blocks[i]->_ivar;
3248 1512 : unsigned int jvar = blocks[i]->_jvar;
3249 :
3250 : // Dirichlet BCs
3251 1512 : std::vector<numeric_index_type> zero_rows;
3252 : PARALLEL_TRY
3253 : {
3254 1512 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3255 49712 : for (const auto & bnode : bnd_nodes)
3256 : {
3257 48200 : BoundaryID boundary_id = bnode->_bnd_id;
3258 48200 : Node * node = bnode->_node;
3259 :
3260 48200 : if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
3261 : {
3262 41544 : const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
3263 :
3264 41544 : if (node->processor_id() == processor_id())
3265 : {
3266 32468 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3267 :
3268 160190 : for (const auto & bc : bcs)
3269 127722 : 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 47568 : zero_rows.push_back(node->dof_number(precond_system.number(), 0, 0));
3274 : }
3275 : }
3276 : }
3277 : }
3278 : }
3279 1512 : PARALLEL_CATCH;
3280 :
3281 1512 : jacobian.close();
3282 :
3283 : // This zeroes the rows corresponding to Dirichlet BCs and puts a 1.0 on the diagonal
3284 1512 : if (ivar == jvar)
3285 1396 : jacobian.zero_rows(zero_rows, 1.0);
3286 : else
3287 116 : jacobian.zero_rows(zero_rows, 0.0);
3288 :
3289 1512 : jacobian.close();
3290 1512 : }
3291 588 : }
3292 :
3293 : void
3294 363355 : NonlinearSystemBase::updateActive(THREAD_ID tid)
3295 : {
3296 363355 : _element_dampers.updateActive(tid);
3297 363355 : _nodal_dampers.updateActive(tid);
3298 363355 : _integrated_bcs.updateActive(tid);
3299 363355 : _dg_kernels.updateActive(tid);
3300 363355 : _interface_kernels.updateActive(tid);
3301 363355 : _dirac_kernels.updateActive(tid);
3302 363355 : _kernels.updateActive(tid);
3303 363355 : _nodal_kernels.updateActive(tid);
3304 363355 : if (tid == 0)
3305 : {
3306 333846 : _general_dampers.updateActive();
3307 333846 : _nodal_bcs.updateActive();
3308 333846 : _preset_nodal_bcs.updateActive();
3309 333846 : _ad_preset_nodal_bcs.updateActive();
3310 333846 : _constraints.updateActive();
3311 333846 : _scalar_kernels.updateActive();
3312 : }
3313 363355 : }
3314 :
3315 : Real
3316 1420 : NonlinearSystemBase::computeDamping(const NumericVector<Number> & solution,
3317 : const NumericVector<Number> & update)
3318 : {
3319 : // Default to no damping
3320 1420 : Real damping = 1.0;
3321 1420 : bool has_active_dampers = false;
3322 :
3323 : try
3324 : {
3325 1420 : 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 1370 : 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 1314 : if (_general_dampers.hasActiveObjects())
3355 : {
3356 : PARALLEL_TRY
3357 : {
3358 404 : TIME_SECTION("computeDamping::general", 3, "Computing General Damping");
3359 :
3360 404 : has_active_dampers = true;
3361 404 : const auto & gdampers = _general_dampers.getActiveObjects();
3362 808 : for (const auto & damper : gdampers)
3363 : {
3364 404 : Real gd_damping = damper->computeDamping(solution, update);
3365 : try
3366 : {
3367 404 : damper->checkMinDamping(gd_damping);
3368 : }
3369 8 : catch (MooseException & e)
3370 : {
3371 8 : _fe_problem.setException(e.what());
3372 8 : }
3373 404 : damping = std::min(gd_damping, damping);
3374 : }
3375 404 : }
3376 404 : 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 1420 : _communicator.min(damping);
3387 :
3388 1420 : if (has_active_dampers && damping < 1.0)
3389 995 : _console << " Damping factor: " << damping << std::endl;
3390 :
3391 1420 : return damping;
3392 : }
3393 :
3394 : void
3395 3794438 : NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
3396 : {
3397 3794438 : _fe_problem.clearDiracInfo();
3398 :
3399 3794438 : std::set<const Elem *> dirac_elements;
3400 :
3401 3794438 : if (_dirac_kernels.hasActiveObjects())
3402 : {
3403 38417 : TIME_SECTION("computeDirac", 3, "Computing DiracKernels");
3404 :
3405 : // TODO: Need a threading fix... but it's complicated!
3406 79882 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3407 : {
3408 41481 : const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
3409 96376 : for (const auto & dkernel : dkernels)
3410 : {
3411 54911 : dkernel->clearPoints();
3412 54911 : dkernel->addPoints();
3413 : }
3414 : }
3415 :
3416 38401 : ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
3417 :
3418 38401 : _fe_problem.getDiracElements(dirac_elements);
3419 :
3420 38401 : DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
3421 : // TODO: Make Dirac work thread!
3422 : // Threads::parallel_reduce(range, cd);
3423 :
3424 38401 : cd(range);
3425 38401 : }
3426 3794422 : }
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 360 : NonlinearSystemBase::residualGhosted()
3439 : {
3440 360 : _need_residual_ghosted = true;
3441 360 : if (!_residual_ghosted)
3442 : {
3443 : // The first time we realize we need a ghosted residual vector,
3444 : // we add it.
3445 288 : _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 288 : if (_Re_time)
3454 : {
3455 50 : const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
3456 50 : _Re_time = &system().add_vector(vector_name, false, GHOSTED);
3457 50 : }
3458 288 : if (_Re_non_time)
3459 : {
3460 288 : const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
3461 288 : _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
3462 288 : }
3463 : }
3464 360 : return *_residual_ghosted;
3465 : }
3466 :
3467 : void
3468 66934 : NonlinearSystemBase::augmentSparsity(SparsityPattern::Graph & sparsity,
3469 : std::vector<dof_id_type> & n_nz,
3470 : std::vector<dof_id_type> & n_oz)
3471 : {
3472 66934 : if (_add_implicit_geometric_coupling_entries_to_jacobian)
3473 : {
3474 901 : _fe_problem.updateGeomSearch();
3475 :
3476 901 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
3477 :
3478 901 : findImplicitGeometricCouplingEntries(_fe_problem.geomSearchData(), graph);
3479 :
3480 901 : if (_fe_problem.getDisplacedProblem())
3481 36 : findImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData(),
3482 : graph);
3483 :
3484 901 : const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
3485 901 : 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 901 : const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
3489 901 : const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
3490 :
3491 1758 : for (const auto & git : graph)
3492 : {
3493 857 : dof_id_type dof = git.first;
3494 857 : dof_id_type local_dof = dof - first_dof_on_proc;
3495 :
3496 857 : if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
3497 126 : continue;
3498 :
3499 731 : const auto & row = git.second;
3500 :
3501 731 : SparsityPattern::Row & sparsity_row = sparsity[local_dof];
3502 :
3503 731 : unsigned int original_row_length = sparsity_row.size();
3504 :
3505 731 : sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
3506 :
3507 1462 : SparsityPattern::sort_row(
3508 731 : sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
3509 :
3510 : // Fix up nonzero arrays
3511 4807 : for (const auto & coupled_dof : row)
3512 : {
3513 4076 : 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 3464 : if (n_nz[local_dof] < n_dofs_on_proc)
3521 3464 : n_nz[local_dof]++;
3522 : }
3523 : }
3524 : }
3525 901 : }
3526 66934 : }
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 14426 : NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
3554 : {
3555 14426 : if (_preconditioner.get() != nullptr)
3556 7 : mooseError("More than one active Preconditioner detected");
3557 :
3558 14419 : _preconditioner = pc;
3559 14419 : }
3560 :
3561 : MoosePreconditioner const *
3562 53655 : NonlinearSystemBase::getPreconditioner() const
3563 : {
3564 53655 : return _preconditioner.get();
3565 : }
3566 :
3567 : void
3568 322 : NonlinearSystemBase::setupDampers()
3569 : {
3570 322 : _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
3571 322 : }
3572 :
3573 : void
3574 1753 : NonlinearSystemBase::reinitIncrementAtQpsForDampers(THREAD_ID /*tid*/,
3575 : const std::set<MooseVariable *> & damped_vars)
3576 : {
3577 3529 : for (const auto & var : damped_vars)
3578 1776 : var->computeIncrementAtQps(*_increment_vec);
3579 1753 : }
3580 :
3581 : void
3582 7441 : NonlinearSystemBase::reinitIncrementAtNodeForDampers(THREAD_ID /*tid*/,
3583 : const std::set<MooseVariable *> & damped_vars)
3584 : {
3585 14882 : for (const auto & var : damped_vars)
3586 7441 : var->computeIncrementAtNode(*_increment_vec);
3587 7441 : }
3588 :
3589 : void
3590 42526 : NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
3591 : {
3592 : // Obtain all blocks and variables covered by all kernels
3593 42526 : std::set<SubdomainID> input_subdomains;
3594 42526 : std::set<std::string> kernel_variables;
3595 :
3596 42526 : bool global_kernels_exist = false;
3597 42526 : global_kernels_exist |= _scalar_kernels.hasActiveObjects();
3598 42526 : global_kernels_exist |= _nodal_kernels.hasActiveObjects();
3599 :
3600 42526 : _kernels.subdomainsCovered(input_subdomains, kernel_variables);
3601 42526 : _dg_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3602 42526 : _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3603 42526 : _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3604 42526 : _constraints.subdomainsCovered(input_subdomains, kernel_variables);
3605 :
3606 42526 : if (_fe_problem.haveFV())
3607 : {
3608 3486 : std::vector<FVElementalKernel *> fv_elemental_kernels;
3609 3486 : _fe_problem.theWarehouse()
3610 6972 : .query()
3611 3486 : .template condition<AttribSystem>("FVElementalKernel")
3612 3486 : .queryInto(fv_elemental_kernels);
3613 :
3614 7454 : for (auto fv_kernel : fv_elemental_kernels)
3615 : {
3616 3968 : if (fv_kernel->blockRestricted())
3617 1972 : for (auto block_id : fv_kernel->blockIDs())
3618 1103 : input_subdomains.insert(block_id);
3619 : else
3620 3099 : global_kernels_exist = true;
3621 3968 : kernel_variables.insert(fv_kernel->variable().name());
3622 :
3623 : // Check for lagrange multiplier
3624 3968 : if (dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel))
3625 134 : kernel_variables.insert(dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel)
3626 268 : ->lambdaVariable()
3627 134 : .name());
3628 : }
3629 :
3630 3486 : std::vector<FVFluxKernel *> fv_flux_kernels;
3631 3486 : _fe_problem.theWarehouse()
3632 6972 : .query()
3633 3486 : .template condition<AttribSystem>("FVFluxKernel")
3634 3486 : .queryInto(fv_flux_kernels);
3635 :
3636 8713 : for (auto fv_kernel : fv_flux_kernels)
3637 : {
3638 5227 : if (fv_kernel->blockRestricted())
3639 2350 : for (auto block_id : fv_kernel->blockIDs())
3640 1250 : input_subdomains.insert(block_id);
3641 : else
3642 4127 : global_kernels_exist = true;
3643 5227 : kernel_variables.insert(fv_kernel->variable().name());
3644 : }
3645 :
3646 3486 : std::vector<FVInterfaceKernel *> fv_interface_kernels;
3647 3486 : _fe_problem.theWarehouse()
3648 6972 : .query()
3649 3486 : .template condition<AttribSystem>("FVInterfaceKernel")
3650 3486 : .queryInto(fv_interface_kernels);
3651 :
3652 3791 : for (auto fvik : fv_interface_kernels)
3653 305 : if (auto scalar_fvik = dynamic_cast<FVScalarLagrangeMultiplierInterface *>(fvik))
3654 14 : kernel_variables.insert(scalar_fvik->lambdaVariable().name());
3655 :
3656 3486 : std::vector<FVFluxBC *> fv_flux_bcs;
3657 3486 : _fe_problem.theWarehouse()
3658 6972 : .query()
3659 3486 : .template condition<AttribSystem>("FVFluxBC")
3660 3486 : .queryInto(fv_flux_bcs);
3661 :
3662 5502 : for (auto fvbc : fv_flux_bcs)
3663 2016 : if (auto scalar_fvbc = dynamic_cast<FVBoundaryScalarLagrangeMultiplierConstraint *>(fvbc))
3664 14 : kernel_variables.insert(scalar_fvbc->lambdaVariable().name());
3665 3486 : }
3666 :
3667 : // Check kernel coverage of subdomains (blocks) in your mesh
3668 42526 : if (!global_kernels_exist)
3669 : {
3670 38575 : std::set<SubdomainID> difference;
3671 38575 : 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 38813 : for (const auto & id : _mesh.interiorLowerDBlocks())
3679 238 : difference.erase(id);
3680 38813 : for (const auto & id : _mesh.boundaryLowerDBlocks())
3681 238 : difference.erase(id);
3682 :
3683 38575 : 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 38567 : }
3705 :
3706 : // Check kernel use of variables
3707 42518 : std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
3708 :
3709 42518 : std::set<VariableName> difference;
3710 42518 : 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 42518 : std::set<VariableName> vars(difference);
3718 42881 : for (auto & var_name : vars)
3719 : {
3720 363 : auto blks = getSubdomainsForVar(var_name);
3721 744 : for (const auto & id : blks)
3722 381 : if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
3723 381 : difference.erase(var_name);
3724 363 : }
3725 :
3726 42518 : 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 42510 : }
3737 :
3738 : bool
3739 28665 : NonlinearSystemBase::containsTimeKernel()
3740 : {
3741 28665 : auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3742 :
3743 28665 : 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 27515 : NonlinearSystemBase::needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3760 : {
3761 27515 : return _integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid);
3762 : }
3763 :
3764 : bool
3765 1282 : NonlinearSystemBase::needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3766 : {
3767 1282 : return _interface_kernels.hasActiveBoundaryObjects(bnd_id, tid);
3768 : }
3769 :
3770 : bool
3771 12193 : NonlinearSystemBase::needSubdomainMaterialOnSide(SubdomainID /*subdomain_id*/,
3772 : THREAD_ID /*tid*/) const
3773 : {
3774 12193 : return _doing_dg;
3775 : }
3776 :
3777 : bool
3778 0 : NonlinearSystemBase::doingDG() const
3779 : {
3780 0 : return _doing_dg;
3781 : }
3782 :
3783 : void
3784 532 : NonlinearSystemBase::setPreviousNewtonSolution(const NumericVector<Number> & soln)
3785 : {
3786 532 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
3787 532 : getVector(Moose::PREVIOUS_NL_SOLUTION_TAG) = soln;
3788 532 : }
3789 :
3790 : void
3791 3798519 : 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 3807551 : for (auto & map_pr : _undisplaced_mortar_functors)
3800 9032 : map_pr.second(compute_type, vector_tags, matrix_tags);
3801 :
3802 3803120 : for (auto & map_pr : _displaced_mortar_functors)
3803 4601 : 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 3798519 : }
3813 :
3814 : void
3815 439 : NonlinearSystemBase::setupScalingData()
3816 : {
3817 439 : 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 439 : const auto n_vars = system().n_vars();
3822 :
3823 439 : if (_scaling_group_variables.empty())
3824 : {
3825 427 : _var_to_group_var.reserve(n_vars);
3826 427 : _num_scaling_groups = n_vars;
3827 :
3828 1128 : for (const auto var_number : make_range(n_vars))
3829 701 : _var_to_group_var.emplace(var_number, var_number);
3830 : }
3831 : else
3832 : {
3833 12 : std::set<unsigned int> var_numbers, var_numbers_covered, var_numbers_not_covered;
3834 48 : for (const auto var_number : make_range(n_vars))
3835 36 : var_numbers.insert(var_number);
3836 :
3837 12 : _num_scaling_groups = _scaling_group_variables.size();
3838 :
3839 24 : for (const auto group_index : index_range(_scaling_group_variables))
3840 48 : for (const auto & var_name : _scaling_group_variables[group_index])
3841 : {
3842 36 : 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 36 : hasVariable(var_name)
3850 24 : ? static_cast<MooseVariableBase &>(getVariable(0, var_name))
3851 60 : : static_cast<MooseVariableBase &>(getScalarVariable(0, var_name));
3852 36 : auto map_pair = _var_to_group_var.emplace(var.number(), group_index);
3853 36 : if (!map_pair.second)
3854 0 : mooseError("Variable ", var_name, " is contained in multiple scaling grouplings");
3855 36 : var_numbers_covered.insert(var.number());
3856 : }
3857 :
3858 12 : 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 12 : _num_scaling_groups = _scaling_group_variables.size() + var_numbers_not_covered.size();
3865 :
3866 12 : auto index = static_cast<unsigned int>(_scaling_group_variables.size());
3867 12 : for (auto var_number : var_numbers_not_covered)
3868 0 : _var_to_group_var.emplace(var_number, index++);
3869 12 : }
3870 :
3871 439 : _variable_autoscaled.resize(n_vars, true);
3872 439 : const auto & number_to_var_map = _vars[0].numberToVariableMap();
3873 :
3874 439 : 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 439 : _auto_scaling_initd = true;
3883 : }
3884 :
3885 : bool
3886 1112 : NonlinearSystemBase::computeScaling()
3887 : {
3888 1112 : if (_compute_scaling_once && _computed_scaling)
3889 460 : return true;
3890 :
3891 652 : _console << "\nPerforming automatic scaling calculation\n" << std::endl;
3892 :
3893 652 : 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 652 : assembleScalingVector();
3898 :
3899 : // container for repeated access of element global dof indices
3900 652 : std::vector<dof_id_type> dof_indices;
3901 :
3902 652 : if (!_auto_scaling_initd)
3903 439 : setupScalingData();
3904 :
3905 652 : std::vector<Real> inverse_scaling_factors(_num_scaling_groups, 0);
3906 652 : std::vector<Real> resid_inverse_scaling_factors(_num_scaling_groups, 0);
3907 652 : std::vector<Real> jac_inverse_scaling_factors(_num_scaling_groups, 0);
3908 652 : auto & dof_map = dofMap();
3909 :
3910 : // what types of scaling do we want?
3911 652 : bool jac_scaling = _resid_vs_jac_scaling_param < 1. - TOLERANCE;
3912 652 : bool resid_scaling = _resid_vs_jac_scaling_param > TOLERANCE;
3913 :
3914 652 : const NumericVector<Number> & scaling_residual = RHS();
3915 :
3916 652 : 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 602 : auto init_vector = NumericVector<Number>::build(this->comm());
3926 602 : init_vector->init(system().n_dofs(), system().n_local_dofs(), /*fast=*/false, PARALLEL);
3927 :
3928 602 : _scaling_matrix->clear();
3929 602 : _scaling_matrix->init(*init_vector);
3930 602 : }
3931 :
3932 602 : _fe_problem.computingScalingJacobian(true);
3933 : // Dispatch to derived classes to ensure that we use the correct matrix tag
3934 602 : computeScalingJacobian();
3935 602 : _fe_problem.computingScalingJacobian(false);
3936 : }
3937 :
3938 652 : 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 652 : if (_fe_problem.getFailNextNonlinearConvergenceCheck())
3950 0 : return false;
3951 :
3952 228330 : 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 4520859 : &scaling_residual](const auto & dof_indices, const auto var_number)
3959 : {
3960 1136309 : for (auto dof_index : dof_indices)
3961 907979 : if (dof_map.local_index(dof_index))
3962 : {
3963 903220 : if (jac_scaling)
3964 : {
3965 : // For now we will use the diagonal for determining scaling
3966 899103 : auto mat_value = (*_scaling_matrix)(dof_index, dof_index);
3967 899103 : auto & factor = jac_inverse_scaling_factors[_var_to_group_var[var_number]];
3968 899103 : factor = std::max(factor, std::abs(mat_value));
3969 : }
3970 903220 : 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 228330 : };
3978 :
3979 : // Compute our scaling factors for the spatial field variables
3980 82245 : for (const auto & elem : _fe_problem.getCurrentAlgebraicElementRange())
3981 310836 : for (const auto i : make_range(system().n_vars()))
3982 229243 : if (_variable_autoscaled[i] && system().variable_type(i).family != SCALAR)
3983 : {
3984 228237 : dof_map.dof_indices(elem, dof_indices, i);
3985 228237 : examine_dof_indices(dof_indices, i);
3986 : }
3987 :
3988 1712 : for (const auto i : make_range(system().n_vars()))
3989 1060 : if (_variable_autoscaled[i] && system().variable_type(i).family == SCALAR)
3990 : {
3991 93 : dof_map.SCALAR_dof_indices(dof_indices, i);
3992 93 : examine_dof_indices(dof_indices, i);
3993 : }
3994 :
3995 652 : if (resid_scaling)
3996 50 : _communicator.max(resid_inverse_scaling_factors);
3997 652 : if (jac_scaling)
3998 602 : _communicator.max(jac_inverse_scaling_factors);
3999 :
4000 652 : 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 652 : else if (jac_scaling)
4020 602 : 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 1688 : for (auto & scaling_factor : inverse_scaling_factors)
4028 1036 : if (scaling_factor == 0)
4029 40 : scaling_factor = 1;
4030 :
4031 : // Now flatten the group scaling factors to the individual variable scaling factors
4032 652 : std::vector<Real> flattened_inverse_scaling_factors(system().n_vars());
4033 1712 : for (const auto i : index_range(flattened_inverse_scaling_factors))
4034 1060 : flattened_inverse_scaling_factors[i] = inverse_scaling_factors[_var_to_group_var[i]];
4035 :
4036 : // Now set the scaling factors for the variables
4037 652 : applyScalingFactors(flattened_inverse_scaling_factors);
4038 652 : if (auto displaced_problem = _fe_problem.getDisplacedProblem().get())
4039 63 : displaced_problem->systemBaseNonlinear(number()).applyScalingFactors(
4040 : flattened_inverse_scaling_factors);
4041 :
4042 652 : _computed_scaling = true;
4043 652 : return true;
4044 652 : }
4045 :
4046 : void
4047 313141 : NonlinearSystemBase::assembleScalingVector()
4048 : {
4049 313141 : if (!hasVector("scaling_factors"))
4050 : // No variables have indicated they need scaling
4051 312591 : return;
4052 :
4053 550 : auto & scaling_vector = getVector("scaling_factors");
4054 :
4055 550 : const auto & lm_mesh = _mesh.getMesh();
4056 550 : const auto & dof_map = dofMap();
4057 :
4058 550 : const auto & field_variables = _vars[0].fieldVariables();
4059 550 : const auto & scalar_variables = _vars[0].scalars();
4060 :
4061 550 : std::vector<dof_id_type> dof_indices;
4062 :
4063 550 : for (const Elem * const elem :
4064 26244 : as_range(lm_mesh.active_local_elements_begin(), lm_mesh.active_local_elements_end()))
4065 31020 : for (const auto * const field_var : field_variables)
4066 : {
4067 18448 : const auto & factors = field_var->arrayScalingFactor();
4068 36896 : for (const auto i : make_range(field_var->count()))
4069 : {
4070 18448 : dof_map.dof_indices(elem, dof_indices, field_var->number() + i);
4071 76448 : for (const auto dof : dof_indices)
4072 58000 : scaling_vector.set(dof, factors[i]);
4073 : }
4074 550 : }
4075 :
4076 662 : 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 112 : dof_map.SCALAR_dof_indices(dof_indices, scalar_var->number());
4081 224 : for (const auto dof : dof_indices)
4082 112 : scaling_vector.set(dof, scalar_var->scalingFactor());
4083 : }
4084 :
4085 : // Parallel assemble
4086 550 : scaling_vector.close();
4087 :
4088 550 : 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 550 : }
4092 :
4093 : bool
4094 312489 : NonlinearSystemBase::preSolve()
4095 : {
4096 : // Clear the iteration counters
4097 312489 : _current_l_its.clear();
4098 312489 : _current_nl_its = 0;
4099 :
4100 : // Initialize the solution vector using a predictor and known values from nodal bcs
4101 312489 : setInitialSolution();
4102 :
4103 : // Now that the initial solution has ben set, potentially perform a residual/Jacobian evaluation
4104 : // to determine variable scaling factors
4105 312489 : if (_automatic_scaling)
4106 : {
4107 1112 : const bool scaling_succeeded = computeScaling();
4108 1112 : 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 312489 : assembleScalingVector();
4115 :
4116 312489 : return true;
4117 : }
4118 :
4119 : void
4120 5612 : NonlinearSystemBase::destroyColoring()
4121 : {
4122 5612 : if (matrixFromColoring())
4123 117 : LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
4124 5612 : }
|
