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