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