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