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