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