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 606320 : 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 606320 : std::vector<T *> objects;
131 : warehouse.query()
132 1212640 : .template condition<AttribSystem>(system_name)
133 606320 : .template condition<AttribSysNum>(system_number)
134 606320 : .template condition<AttribThread>(tid)
135 606320 : .queryInto(objects);
136 :
137 1164451 : for (auto * object : objects)
138 558131 : results.push_back(object);
139 606320 : }
140 : }
141 :
142 62361 : NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem,
143 : System & sys,
144 62361 : const std::string & name)
145 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
146 62361 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
147 62361 : _sys(sys),
148 62361 : _last_nl_rnorm(0.),
149 62361 : _current_nl_its(0),
150 62361 : _residual_ghosted(NULL),
151 62361 : _Re_time_tag(-1),
152 62361 : _Re_time(NULL),
153 62361 : _Re_non_time_tag(-1),
154 62361 : _Re_non_time(NULL),
155 62361 : _scalar_kernels(/*threaded=*/false),
156 62361 : _nodal_bcs(/*threaded=*/false),
157 62361 : _preset_nodal_bcs(/*threaded=*/false),
158 62361 : _ad_preset_nodal_bcs(/*threaded=*/false),
159 : #ifdef MOOSE_KOKKOS_ENABLED
160 47279 : _kokkos_kernels(/*threaded=*/false),
161 47279 : _kokkos_integrated_bcs(/*threaded=*/false),
162 47279 : _kokkos_nodal_bcs(/*threaded=*/false),
163 47279 : _kokkos_preset_nodal_bcs(/*threaded=*/false),
164 47279 : _kokkos_nodal_kernels(/*threaded=*/false),
165 : #endif
166 62361 : _general_dampers(/*threaded=*/false),
167 62361 : _splits(/*threaded=*/false),
168 62361 : _increment_vec(NULL),
169 62361 : _use_finite_differenced_preconditioner(false),
170 62361 : _fdcoloring(nullptr),
171 62361 : _fsp(nullptr),
172 62361 : _add_implicit_geometric_coupling_entries_to_jacobian(false),
173 62361 : _assemble_constraints_separately(false),
174 62361 : _need_residual_ghosted(false),
175 62361 : _debugging_residuals(false),
176 62361 : _doing_dg(false),
177 62361 : _n_iters(0),
178 62361 : _n_linear_iters(0),
179 62361 : _n_residual_evaluations(0),
180 62361 : _final_residual(0.),
181 62361 : _computing_pre_smo_residual(false),
182 62361 : _pre_smo_residual(0),
183 62361 : _initial_residual(0),
184 62361 : _use_pre_smo_residual(false),
185 62361 : _print_all_var_norms(false),
186 62361 : _has_save_in(false),
187 62361 : _has_diag_save_in(false),
188 62361 : _has_nodalbc_save_in(false),
189 62361 : _has_nodalbc_diag_save_in(false),
190 62361 : _computed_scaling(false),
191 62361 : _compute_scaling_once(true),
192 62361 : _resid_vs_jac_scaling_param(0),
193 62361 : _off_diagonals_in_auto_scaling(false),
194 498888 : _auto_scaling_initd(false)
195 : {
196 62361 : getResidualNonTimeVector();
197 : // Don't need to add the matrix - it already exists (for now)
198 62361 : _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 62361 : _fe_problem.addMatrixTag("TIME");
203 :
204 62361 : _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
205 :
206 62361 : _sys.identify_variable_groups(_fe_problem.identifyVariableGroupsInNL());
207 :
208 62361 : if (!_fe_problem.defaultGhosting())
209 : {
210 62288 : auto & dof_map = _sys.get_dof_map();
211 62288 : dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting());
212 62288 : dof_map.set_implicit_neighbor_dofs(false);
213 : }
214 62361 : }
215 :
216 59201 : NonlinearSystemBase::~NonlinearSystemBase() = default;
217 :
218 : void
219 60500 : NonlinearSystemBase::preInit()
220 : {
221 60500 : SolverSystem::preInit();
222 :
223 60500 : if (_fe_problem.hasDampers())
224 159 : setupDampers();
225 :
226 60500 : if (_residual_copy.get())
227 0 : _residual_copy->init(_sys.n_dofs(), false, SERIAL);
228 :
229 : #ifdef MOOSE_KOKKOS_ENABLED
230 45870 : if (_fe_problem.hasKokkosResidualObjects())
231 2039 : _sys.get_dof_map().full_sparsity_pattern_needed();
232 : #endif
233 60500 : }
234 :
235 : void
236 7016 : NonlinearSystemBase::reinitMortarFunctors()
237 : {
238 : // reinit is called on meshChanged() in FEProblemBase. We could implement meshChanged() instead.
239 : // Subdomains might have changed
240 7086 : for (auto & functor : _displaced_mortar_functors)
241 70 : functor.second.setupMortarMaterials();
242 7118 : for (auto & functor : _undisplaced_mortar_functors)
243 102 : functor.second.setupMortarMaterials();
244 7016 : }
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 121264 : NonlinearSystemBase::getFVSetupObjects(THREAD_ID tid)
255 : {
256 121264 : std::vector<SetupInterface *> fv_objects;
257 121264 : auto & warehouse = _fe_problem.theWarehouse();
258 :
259 242528 : appendFVSetupObjects<FVElementalKernel>(
260 : warehouse, "FVElementalKernel", number(), tid, fv_objects);
261 242528 : appendFVSetupObjects<FVFluxKernel>(warehouse, "FVFluxKernel", number(), tid, fv_objects);
262 242528 : appendFVSetupObjects<FVBoundaryCondition>(warehouse, "FVDirichletBC", number(), tid, fv_objects);
263 242528 : appendFVSetupObjects<FVBoundaryCondition>(warehouse, "FVFluxBC", number(), tid, fv_objects);
264 242528 : appendFVSetupObjects<FVInterfaceKernel>(
265 : warehouse, "FVInterfaceKernel", number(), tid, fv_objects);
266 :
267 121264 : return fv_objects;
268 0 : }
269 :
270 : void
271 59766 : NonlinearSystemBase::initialSetup()
272 : {
273 298830 : TIME_SECTION("nlInitialSetup", 2, "Setting Up Nonlinear System");
274 :
275 59766 : SolverSystem::initialSetup();
276 :
277 : {
278 298830 : TIME_SECTION("kernelsInitialSetup", 2, "Setting Up Kernels/BCs/Constraints");
279 :
280 125543 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
281 : {
282 65777 : _kernels.initialSetup(tid);
283 65777 : _nodal_kernels.initialSetup(tid);
284 65777 : _dirac_kernels.initialSetup(tid);
285 65777 : if (_doing_dg)
286 1084 : _dg_kernels.initialSetup(tid);
287 65777 : _interface_kernels.initialSetup(tid);
288 :
289 65777 : _element_dampers.initialSetup(tid);
290 65777 : _nodal_dampers.initialSetup(tid);
291 65777 : _integrated_bcs.initialSetup(tid);
292 :
293 65777 : if (_fe_problem.haveFV())
294 16848 : for (auto * fv_object : getFVSetupObjects(tid))
295 16848 : fv_object->initialSetup();
296 : }
297 :
298 59766 : _scalar_kernels.initialSetup();
299 59766 : _constraints.initialSetup();
300 59766 : _general_dampers.initialSetup();
301 59766 : _nodal_bcs.initialSetup();
302 59766 : _preset_nodal_bcs.residualSetup();
303 59766 : _ad_preset_nodal_bcs.residualSetup();
304 :
305 : #ifdef MOOSE_KOKKOS_ENABLED
306 45369 : _kokkos_kernels.initialSetup();
307 45369 : _kokkos_nodal_kernels.initialSetup();
308 45369 : _kokkos_integrated_bcs.initialSetup();
309 45369 : _kokkos_nodal_bcs.initialSetup();
310 : #endif
311 59766 : }
312 :
313 : {
314 298830 : TIME_SECTION("mortarSetup", 2, "Initializing Mortar Interfaces");
315 :
316 119532 : auto create_mortar_functors = [this](const bool displaced)
317 : {
318 : // go over mortar interfaces and construct functors
319 119532 : const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
320 120596 : 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 119532 : };
344 :
345 59766 : create_mortar_functors(false);
346 59766 : create_mortar_functors(true);
347 59766 : }
348 :
349 59766 : 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 59766 : if (_preconditioner)
358 13288 : _preconditioner->initialSetup();
359 59766 : }
360 :
361 : void
362 273114 : NonlinearSystemBase::timestepSetup()
363 : {
364 273114 : SolverSystem::timestepSetup();
365 :
366 573364 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
367 : {
368 300250 : _kernels.timestepSetup(tid);
369 300250 : _nodal_kernels.timestepSetup(tid);
370 300250 : _dirac_kernels.timestepSetup(tid);
371 300250 : if (_doing_dg)
372 1609 : _dg_kernels.timestepSetup(tid);
373 300250 : _interface_kernels.timestepSetup(tid);
374 300250 : _element_dampers.timestepSetup(tid);
375 300250 : _nodal_dampers.timestepSetup(tid);
376 300250 : _integrated_bcs.timestepSetup(tid);
377 :
378 300250 : if (_fe_problem.haveFV())
379 80379 : for (auto * fv_object : getFVSetupObjects(tid))
380 80379 : fv_object->timestepSetup();
381 : }
382 273114 : _scalar_kernels.timestepSetup();
383 273114 : _constraints.timestepSetup();
384 273114 : _general_dampers.timestepSetup();
385 273114 : _nodal_bcs.timestepSetup();
386 273114 : _preset_nodal_bcs.timestepSetup();
387 273114 : _ad_preset_nodal_bcs.timestepSetup();
388 :
389 : #ifdef MOOSE_KOKKOS_ENABLED
390 200560 : _kokkos_kernels.timestepSetup();
391 200560 : _kokkos_nodal_kernels.timestepSetup();
392 200560 : _kokkos_integrated_bcs.timestepSetup();
393 200560 : _kokkos_nodal_bcs.timestepSetup();
394 : #endif
395 273114 : }
396 :
397 : void
398 1808570 : NonlinearSystemBase::customSetup(const ExecFlagType & exec_type)
399 : {
400 1808570 : SolverSystem::customSetup(exec_type);
401 :
402 3799245 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
403 : {
404 1990675 : _kernels.customSetup(exec_type, tid);
405 1990675 : _nodal_kernels.customSetup(exec_type, tid);
406 1990675 : _dirac_kernels.customSetup(exec_type, tid);
407 1990675 : if (_doing_dg)
408 11026 : _dg_kernels.customSetup(exec_type, tid);
409 1990675 : _interface_kernels.customSetup(exec_type, tid);
410 1990675 : _element_dampers.customSetup(exec_type, tid);
411 1990675 : _nodal_dampers.customSetup(exec_type, tid);
412 1990675 : _integrated_bcs.customSetup(exec_type, tid);
413 :
414 1990675 : if (_fe_problem.haveFV())
415 582168 : for (auto * fv_object : getFVSetupObjects(tid))
416 582168 : fv_object->customSetup(exec_type);
417 : }
418 1808570 : _scalar_kernels.customSetup(exec_type);
419 1808570 : _constraints.customSetup(exec_type);
420 1808570 : _general_dampers.customSetup(exec_type);
421 1808570 : _nodal_bcs.customSetup(exec_type);
422 1808570 : _preset_nodal_bcs.customSetup(exec_type);
423 1808570 : _ad_preset_nodal_bcs.customSetup(exec_type);
424 :
425 : #ifdef MOOSE_KOKKOS_ENABLED
426 1323856 : _kokkos_kernels.customSetup(exec_type);
427 1323856 : _kokkos_nodal_kernels.customSetup(exec_type);
428 1323856 : _kokkos_integrated_bcs.customSetup(exec_type);
429 1323856 : _kokkos_nodal_bcs.customSetup(exec_type);
430 : #endif
431 1808570 : }
432 :
433 : void
434 323334 : NonlinearSystemBase::setupDM()
435 : {
436 323334 : if (_fsp)
437 113 : _fsp->setupDM();
438 323334 : }
439 :
440 : void
441 77437 : NonlinearSystemBase::addKernel(const std::string & kernel_name,
442 : const std::string & name,
443 : InputParameters & parameters)
444 : {
445 162615 : 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 85331 : _factory.create<KernelBase>(kernel_name, name, parameters, tid);
450 85196 : _kernels.addObject(kernel, tid);
451 85181 : postAddResidualObject(*kernel);
452 : // Add to theWarehouse, a centralized storage for all moose objects
453 85178 : _fe_problem.theWarehouse().add(kernel);
454 85178 : }
455 :
456 77284 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
457 120 : _has_save_in = true;
458 77284 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
459 96 : _has_diag_save_in = true;
460 77284 : }
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 75520 : NonlinearSystemBase::addBoundaryCondition(const std::string & bc_name,
518 : const std::string & name,
519 : InputParameters & parameters)
520 : {
521 : // ThreadID
522 75520 : THREAD_ID tid = 0;
523 :
524 : // Create the object
525 : std::shared_ptr<BoundaryCondition> bc =
526 75520 : _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
527 75479 : postAddResidualObject(*bc);
528 :
529 : // Active BoundaryIDs for the object
530 75479 : const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
531 75479 : auto bc_var = dynamic_cast<const MooseVariableFieldBase *>(&bc->variable());
532 75479 : _vars[tid].addBoundaryVar(boundary_ids, bc_var);
533 :
534 : // Cast to the various types of BCs
535 75479 : std::shared_ptr<NodalBCBase> nbc = std::dynamic_pointer_cast<NodalBCBase>(bc);
536 75479 : std::shared_ptr<IntegratedBCBase> ibc = std::dynamic_pointer_cast<IntegratedBCBase>(bc);
537 :
538 : // NodalBCBase
539 75479 : if (nbc)
540 : {
541 67161 : 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 67158 : _nodal_bcs.addObject(nbc);
549 : // Add to theWarehouse, a centralized storage for all moose objects
550 67158 : _fe_problem.theWarehouse().add(nbc);
551 67158 : _vars[tid].addBoundaryVars(boundary_ids, nbc->getCoupledVars());
552 :
553 67158 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
554 90 : _has_nodalbc_save_in = true;
555 67158 : 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 67158 : std::shared_ptr<DirichletBCBase> dbc = std::dynamic_pointer_cast<DirichletBCBase>(bc);
560 67158 : if (dbc && dbc->preset())
561 57749 : _preset_nodal_bcs.addObject(dbc);
562 :
563 67158 : std::shared_ptr<ADDirichletBCBase> addbc = std::dynamic_pointer_cast<ADDirichletBCBase>(bc);
564 67158 : if (addbc && addbc->preset())
565 1593 : _ad_preset_nodal_bcs.addObject(addbc);
566 67158 : }
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 75476 : }
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 289811 : NonlinearSystemBase::shouldEvaluatePreSMOResidual() const
725 : {
726 289811 : 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 279876 : if (_app.parameters().get<bool>("use_legacy_initial_residual_evaluation_behavior"))
736 16 : return true;
737 :
738 279860 : return _use_pre_smo_residual;
739 : }
740 :
741 : Real
742 752238 : NonlinearSystemBase::referenceResidual() const
743 : {
744 752238 : 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 752249 : NonlinearSystemBase::initialResidual() const
758 : {
759 752249 : return _initial_residual;
760 : }
761 :
762 : void
763 286924 : NonlinearSystemBase::setInitialResidual(Real r)
764 : {
765 286924 : _initial_residual = r;
766 286924 : }
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 3055432 : NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
802 : {
803 : parallel_object_only();
804 :
805 9166296 : TIME_SECTION("nl::computeResidualTags", 5);
806 :
807 3055432 : _fe_problem.setCurrentNonlinearSystem(number());
808 3055432 : _fe_problem.setCurrentlyComputingResidual(true);
809 :
810 3055432 : bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
811 :
812 3055432 : _n_residual_evaluations++;
813 :
814 : // not suppose to do anythin on matrix
815 3055432 : deactivateAllMatrixTags();
816 :
817 3055432 : FloatingPointExceptionGuard fpe_guard(_app);
818 :
819 3055432 : 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 3055432 : zeroTaggedVectors(tags);
830 3055432 : computeResidualInternal(tags);
831 3055293 : closeTaggedVectors(tags);
832 :
833 3055293 : if (required_residual)
834 : {
835 3024598 : auto & residual = getVector(residualVectorTag());
836 3024598 : if (!_time_integrators.empty())
837 : {
838 5281193 : for (auto & ti : _time_integrators)
839 2647712 : ti->postResidual(residual);
840 : }
841 : else
842 391117 : residual += *_Re_non_time;
843 3024598 : residual.close();
844 : }
845 3055293 : if (_fe_problem.computingScalingResidual())
846 : // We don't want to do nodal bcs or anything else
847 45 : return;
848 :
849 3055248 : computeNodalBCsResidual(tags);
850 3055248 : closeTaggedVectors(tags);
851 :
852 : // If we are debugging residuals we need one more assignment to have the ghosted copy up to
853 : // date
854 3055248 : 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 3055248 : if (_has_nodalbc_save_in)
863 157 : _fe_problem.getAuxiliarySystem().solution().close();
864 3055248 : 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 3055357 : activateAllMatrixTags();
876 :
877 3055357 : _fe_problem.setCurrentlyComputingResidual(false);
878 3055447 : }
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 243843 : NonlinearSystemBase::onTimestepBegin()
921 : {
922 488026 : for (auto & ti : _time_integrators)
923 244183 : ti->preSolve();
924 243843 : if (_predictor.get())
925 353 : _predictor->timestepSetup();
926 243843 : }
927 :
928 : void
929 290807 : NonlinearSystemBase::setInitialSolution()
930 : {
931 290807 : deactivateAllMatrixTags();
932 :
933 290807 : NumericVector<Number> & initial_solution(solution());
934 290807 : 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 1454035 : TIME_SECTION("initialBCs", 2, "Applying BCs To Initial Condition");
950 :
951 290807 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
952 14411366 : for (const auto & bnode : bnd_nodes)
953 : {
954 14120559 : BoundaryID boundary_id = bnode->_bnd_id;
955 14120559 : Node * node = bnode->_node;
956 :
957 14120559 : if (node->processor_id() == processor_id())
958 : {
959 10481907 : bool has_preset_nodal_bcs = _preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
960 10481907 : bool has_ad_preset_nodal_bcs = _ad_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
961 :
962 : // reinit variables in nodes
963 10481907 : if (has_preset_nodal_bcs || has_ad_preset_nodal_bcs)
964 3252529 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
965 :
966 10481907 : 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 10481907 : 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 290807 : }
981 :
982 : #ifdef MOOSE_KOKKOS_ENABLED
983 212225 : if (_kokkos_preset_nodal_bcs.hasObjects())
984 4598 : setKokkosInitialSolution();
985 : #endif
986 :
987 290807 : _sys.solution->close();
988 290807 : update();
989 :
990 : // Set constraint secondary values
991 290807 : setConstraintSecondaryValues(initial_solution, false);
992 :
993 290807 : if (_fe_problem.getDisplacedProblem())
994 32738 : setConstraintSecondaryValues(initial_solution, true);
995 290807 : }
996 :
997 : void
998 111 : NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
999 : {
1000 111 : _predictor = predictor;
1001 111 : }
1002 :
1003 : void
1004 5235034 : NonlinearSystemBase::subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
1005 : {
1006 5235034 : SolverSystem::subdomainSetup();
1007 :
1008 5235034 : _kernels.subdomainSetup(subdomain, tid);
1009 5235034 : _nodal_kernels.subdomainSetup(subdomain, tid);
1010 5235034 : _element_dampers.subdomainSetup(subdomain, tid);
1011 5235034 : _nodal_dampers.subdomainSetup(subdomain, tid);
1012 5235034 : }
1013 :
1014 : NumericVector<Number> &
1015 30715 : NonlinearSystemBase::getResidualTimeVector()
1016 : {
1017 30715 : if (!_Re_time)
1018 : {
1019 30637 : _Re_time_tag = _fe_problem.addVectorTag("TIME");
1020 :
1021 : // Most applications don't need the expense of ghosting
1022 30637 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1023 30637 : _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 30715 : return *_Re_time;
1035 : }
1036 :
1037 : NumericVector<Number> &
1038 93266 : NonlinearSystemBase::getResidualNonTimeVector()
1039 : {
1040 93266 : if (!_Re_non_time)
1041 : {
1042 62361 : _Re_non_time_tag = _fe_problem.addVectorTag("NONTIME");
1043 :
1044 : // Most applications don't need the expense of ghosting
1045 62361 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1046 62361 : _Re_non_time = &addVector(_Re_non_time_tag, false, ptype);
1047 : }
1048 30905 : 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 93266 : 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 21962 : NonlinearSystemBase::enforceNodalConstraintsResidual(NumericVector<Number> & residual)
1079 : {
1080 21962 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1081 21962 : residual.close();
1082 21962 : if (_constraints.hasActiveNodalConstraints())
1083 : {
1084 858 : const auto & ncs = _constraints.getActiveNodalConstraints();
1085 1716 : for (const auto & nc : ncs)
1086 : {
1087 858 : std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1088 858 : std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1089 :
1090 858 : if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1091 : {
1092 854 : _fe_problem.reinitNodes(primary_node_ids, tid);
1093 854 : _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1094 854 : nc->computeResidual(residual);
1095 : }
1096 : }
1097 858 : _fe_problem.addCachedResidualDirectly(residual, tid);
1098 858 : residual.close();
1099 : }
1100 21962 : }
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 323545 : 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 356283 : : static_cast<SubProblem &>(_fe_problem);
1205 323545 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1206 :
1207 323545 : bool constraints_applied = false;
1208 :
1209 378763 : 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 323545 : std::set<dof_id_type> unique_secondary_node_ids;
1279 :
1280 727006 : for (const auto & secondary_id : _mesh.meshSubdomains())
1281 : {
1282 1076586 : for (const auto & primary_id : _mesh.meshSubdomains())
1283 : {
1284 673125 : 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 323545 : _communicator.max(constraints_applied);
1329 :
1330 323545 : if (constraints_applied)
1331 : {
1332 692 : solution.close();
1333 692 : update();
1334 : }
1335 323545 : }
1336 :
1337 : void
1338 26409 : NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool displaced)
1339 : {
1340 : // Make sure the residual is in a good state
1341 26409 : 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 30856 : : static_cast<SubProblem &>(_fe_problem);
1349 26409 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1350 :
1351 : bool constraints_applied;
1352 26409 : bool residual_has_inserted_values = false;
1353 26409 : if (!_assemble_constraints_separately)
1354 26409 : constraints_applied = false;
1355 48696 : 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 26409 : if (!_assemble_constraints_separately)
1484 : {
1485 26409 : _communicator.max(constraints_applied);
1486 :
1487 26409 : 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 26409 : THREAD_ID tid = 0;
1505 26409 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
1506 26409 : 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 26409 : std::set<dof_id_type> unique_secondary_node_ids;
1552 :
1553 26409 : constraints_applied = false;
1554 26409 : residual_has_inserted_values = false;
1555 26409 : bool has_writable_variables = false;
1556 90257 : for (const auto & secondary_id : _mesh.meshSubdomains())
1557 : {
1558 272678 : for (const auto & primary_id : _mesh.meshSubdomains())
1559 : {
1560 208830 : 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 26409 : _communicator.max(constraints_applied);
1622 :
1623 26409 : 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 26409 : _communicator.max(has_writable_variables);
1638 :
1639 26409 : 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 26409 : _fe_problem.addCachedResidual(0);
1651 26409 : }
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 3155133 : NonlinearSystemBase::residualSetup()
1704 : {
1705 9465399 : TIME_SECTION("residualSetup", 3);
1706 :
1707 3155133 : SolverSystem::residualSetup();
1708 :
1709 6629176 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1710 : {
1711 3474043 : _kernels.residualSetup(tid);
1712 3474043 : _nodal_kernels.residualSetup(tid);
1713 3474043 : _dirac_kernels.residualSetup(tid);
1714 3474043 : if (_doing_dg)
1715 34549 : _dg_kernels.residualSetup(tid);
1716 3474043 : _interface_kernels.residualSetup(tid);
1717 3474043 : _element_dampers.residualSetup(tid);
1718 3474043 : _nodal_dampers.residualSetup(tid);
1719 3474043 : _integrated_bcs.residualSetup(tid);
1720 : }
1721 3155133 : _scalar_kernels.residualSetup();
1722 3155133 : _constraints.residualSetup();
1723 3155133 : _general_dampers.residualSetup();
1724 3155133 : _nodal_bcs.residualSetup();
1725 3155133 : _preset_nodal_bcs.residualSetup();
1726 3155133 : _ad_preset_nodal_bcs.residualSetup();
1727 :
1728 : #ifdef MOOSE_KOKKOS_ENABLED
1729 2296463 : _kokkos_kernels.residualSetup();
1730 2296463 : _kokkos_nodal_kernels.residualSetup();
1731 2296463 : _kokkos_integrated_bcs.residualSetup();
1732 2296463 : _kokkos_nodal_bcs.residualSetup();
1733 : #endif
1734 :
1735 : // Avoid recursion
1736 3155133 : if (this == &_fe_problem.currentNonlinearSystem())
1737 3065331 : _fe_problem.residualSetup();
1738 3155133 : }
1739 :
1740 : void
1741 3055432 : NonlinearSystemBase::computeResidualInternal(const std::set<TagID> & tags)
1742 : {
1743 : parallel_object_only();
1744 :
1745 9166296 : TIME_SECTION("computeResidualInternal", 3);
1746 :
1747 3055432 : 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 3055432 : std::vector<GeneralUserObject *> uos;
1752 3055432 : _fe_problem.theWarehouse()
1753 6110864 : .query()
1754 3055432 : .condition<AttribSystem>("UserObject")
1755 3055432 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1756 3055432 : .queryInto(uos);
1757 3055432 : for (auto & uo : uos)
1758 0 : uo->residualSetup();
1759 3055432 : for (auto & uo : uos)
1760 : {
1761 0 : uo->initialize();
1762 0 : uo->execute();
1763 0 : uo->finalize();
1764 : }
1765 :
1766 : // reinit scalar variables
1767 6420148 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1768 3364716 : _fe_problem.reinitScalars(tid);
1769 :
1770 : #ifdef MOOSE_KOKKOS_ENABLED
1771 2222896 : if (_fe_problem.hasKokkosResidualObjects())
1772 67208 : computeKokkosResidual(tags);
1773 : #endif
1774 :
1775 : // residual contributions from the domain
1776 : PARALLEL_TRY
1777 : {
1778 9166296 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1779 :
1780 3055432 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
1781 :
1782 3055432 : ComputeResidualThread cr(_fe_problem, tags);
1783 3055432 : 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 3055414 : if (_fe_problem.haveFV())
1791 : {
1792 : ComputeFVFluxResidualThread<FVRange> fvr(
1793 70882 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
1794 70882 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
1795 70882 : Threads::parallel_reduce(faces, fvr);
1796 70882 : }
1797 3055408 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
1798 3055408 : 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 3055408 : }
1806 :
1807 3055408 : unsigned int n_threads = libMesh::n_threads();
1808 6420092 : for (unsigned int i = 0; i < n_threads;
1809 : i++) // Add any cached residuals that might be hanging around
1810 3364684 : _fe_problem.addCachedResidual(i);
1811 3055408 : }
1812 3055408 : 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 3055299 : if (_scalar_kernels.hasActiveObjects())
1819 : {
1820 172047 : 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 57349 : if (!tags.size() || tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL))
1827 57332 : 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 57349 : bool have_scalar_contributions = false;
1836 57349 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
1837 216992 : for (const auto & scalar_kernel : scalars)
1838 : {
1839 159643 : scalar_kernel->reinit();
1840 159643 : const std::vector<dof_id_type> & dof_indices = scalar_kernel->variable().dofIndices();
1841 159643 : const DofMap & dof_map = scalar_kernel->variable().dofMap();
1842 159643 : const dof_id_type first_dof = dof_map.first_dof();
1843 159643 : const dof_id_type end_dof = dof_map.end_dof();
1844 189792 : for (dof_id_type dof : dof_indices)
1845 : {
1846 160220 : if (dof >= first_dof && dof < end_dof)
1847 : {
1848 130071 : scalar_kernel->computeResidual();
1849 130071 : have_scalar_contributions = true;
1850 130071 : break;
1851 : }
1852 : }
1853 : }
1854 57349 : if (have_scalar_contributions)
1855 48426 : _fe_problem.addResidualScalar();
1856 57349 : }
1857 : }
1858 3055299 : PARALLEL_CATCH;
1859 :
1860 : // residual contributions from Block NodalKernels
1861 : PARALLEL_TRY
1862 : {
1863 3055299 : 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 3055299 : PARALLEL_CATCH;
1885 :
1886 3055299 : 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 3055254 : 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 3055254 : PARALLEL_CATCH;
1912 :
1913 3055254 : mortarConstraints(Moose::ComputeType::Residual, tags, {});
1914 :
1915 3055254 : if (_residual_copy.get())
1916 : {
1917 0 : _Re_non_time->close();
1918 0 : _Re_non_time->localize(*_residual_copy);
1919 : }
1920 :
1921 3055254 : 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 3055254 : PARALLEL_TRY { computeDiracContributions(tags, false); }
1929 3055248 : PARALLEL_CATCH;
1930 :
1931 3055248 : if (_fe_problem._has_constraints)
1932 : {
1933 21962 : PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1934 21962 : PARALLEL_CATCH;
1935 21962 : _Re_non_time->close();
1936 : }
1937 :
1938 : // Add in Residual contributions from other Constraints
1939 3055248 : if (_fe_problem._has_constraints)
1940 : {
1941 : PARALLEL_TRY
1942 : {
1943 : // Undisplaced Constraints
1944 21962 : constraintResiduals(*_Re_non_time, false);
1945 :
1946 : // Displaced Constraints
1947 21962 : if (_fe_problem.getDisplacedProblem())
1948 4447 : constraintResiduals(*_Re_non_time, true);
1949 :
1950 21962 : if (_fe_problem.computingNonlinearResid())
1951 9935 : _constraints.residualEnd();
1952 : }
1953 21962 : PARALLEL_CATCH;
1954 21962 : _Re_non_time->close();
1955 : }
1956 :
1957 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
1958 : // counters
1959 3055248 : _app.solutionInvalidity().syncIteration();
1960 3055248 : _app.solutionInvalidity().accumulateIterationIntoTimeStepOccurences();
1961 3055556 : }
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 3055248 : NonlinearSystemBase::computeNodalBCsResidual(const std::set<TagID> & tags)
2090 : {
2091 : #ifdef MOOSE_KOKKOS_ENABLED
2092 2222771 : if (_fe_problem.hasKokkosResidualObjects())
2093 67208 : 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 3055248 : if (_has_save_in)
2099 284 : _fe_problem.getAuxiliarySystem().solution().close();
2100 :
2101 : // Select nodal kernels
2102 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2103 :
2104 3055248 : if (tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL) || !tags.size())
2105 3018833 : nbc_warehouse = &_nodal_bcs;
2106 36415 : else if (tags.size() == 1)
2107 17680 : nbc_warehouse = &(_nodal_bcs.getVectorTagObjectWarehouse(*(tags.begin()), 0));
2108 : else
2109 18735 : nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0));
2110 :
2111 : // Return early if there is no nodal kernel
2112 3055248 : if (!nbc_warehouse->hasActiveObjects())
2113 371634 : return;
2114 :
2115 : PARALLEL_TRY
2116 : {
2117 2683614 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2118 :
2119 2683614 : if (!bnd_nodes.empty())
2120 : {
2121 8049891 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2122 :
2123 140390371 : for (const auto & bnode : bnd_nodes)
2124 : {
2125 137707074 : BoundaryID boundary_id = bnode->_bnd_id;
2126 137707074 : Node * node = bnode->_node;
2127 :
2128 240773114 : if (node->processor_id() == processor_id() &&
2129 103066040 : nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
2130 : {
2131 : // reinit variables in nodes
2132 53195247 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2133 :
2134 53195247 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2135 111343247 : for (const auto & nbc : bcs)
2136 58148000 : if (nbc->shouldApply())
2137 58146773 : nbc->computeResidual();
2138 : }
2139 : }
2140 2683297 : }
2141 : }
2142 2683614 : PARALLEL_CATCH;
2143 :
2144 2683614 : if (_Re_time)
2145 2373799 : _Re_time->close();
2146 2683614 : _Re_non_time->close();
2147 : }
2148 :
2149 : void
2150 475640 : 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 475640 : if (_has_diag_save_in)
2155 170 : _fe_problem.getAuxiliarySystem().solution().close();
2156 :
2157 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2158 :
2159 : // Select nodal kernels
2160 475640 : if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
2161 466401 : 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 475640 : if (!nbc_warehouse->hasActiveObjects())
2169 79563 : 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 1982612 : 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 1586535 : if (nbc_warehouse->hasActiveBoundaryObjects(bid))
2190 : {
2191 777681 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
2192 1633686 : for (const auto & bc : bcs)
2193 : {
2194 856005 : 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 856005 : std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
2199 856927 : 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 856005 : 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 831216 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2212 435139 : _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 20543190 : for (const auto & bnode : bnd_nodes)
2224 : {
2225 20147113 : BoundaryID boundary_id = bnode->_bnd_id;
2226 20147113 : Node * node = bnode->_node;
2227 :
2228 29556713 : if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
2229 9409600 : node->processor_id() == processor_id())
2230 : {
2231 7143119 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2232 :
2233 7143119 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2234 15182540 : for (const auto & bc : bcs)
2235 : {
2236 : // Get the set of involved MOOSE vars for this BC
2237 8039421 : 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 20143411 : for (const auto & it : coupling_entries)
2244 : {
2245 12103990 : 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 12103990 : if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
2252 8043687 : bc->computeOffDiagJacobian(jvar);
2253 : }
2254 :
2255 8039421 : const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars();
2256 8039751 : 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 : // Begin the diagonal node-face constraint accumulation phase for neighbor Jacobian
2546 : // blocks.
2547 3815 : _fe_problem.prepareAssemblyNeighbor(0);
2548 :
2549 3815 : nfc->prepareShapes(nfc->variable().number());
2550 3815 : nfc->prepareNeighborShapes(nfc->variable().number());
2551 :
2552 3815 : nfc->computeJacobian();
2553 :
2554 3815 : if (nfc->overwriteSecondaryJacobian())
2555 : {
2556 : // Add this variable's dof's row to be zeroed
2557 3815 : zero_rows.push_back(nfc->variable().nodalDofIndex());
2558 : }
2559 :
2560 3815 : std::vector<dof_id_type> secondary_dofs(1, nfc->variable().nodalDofIndex());
2561 :
2562 : // Assume that if the user is overwriting the secondary Jacobian, then they are
2563 : // supplying Jacobians that do not correspond to their other physics
2564 : // (e.g. Kernels), hence we should not apply a scalingFactor that is normally
2565 : // based on the order of their other physics (e.g. Kernels)
2566 : Real scaling_factor =
2567 3815 : nfc->overwriteSecondaryJacobian() ? 1. : nfc->variable().scalingFactor();
2568 :
2569 : // Cache the jacobian block for the secondary side
2570 7630 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2571 3815 : nfc->_Kee,
2572 : secondary_dofs,
2573 3815 : nfc->_connected_dof_indices,
2574 : scaling_factor);
2575 :
2576 : // Cache Ken, Kne, Knn
2577 3815 : if (nfc->addCouplingEntriesToJacobian())
2578 : {
2579 : // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2580 : // factor when we're overwriting secondary stuff)
2581 7630 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2582 3815 : nfc->_Ken,
2583 : secondary_dofs,
2584 3815 : nfc->primaryVariable().dofIndicesNeighbor(),
2585 : scaling_factor);
2586 :
2587 : // Use _connected_dof_indices to get all the correct columns
2588 7630 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2589 3815 : nfc->_Kne,
2590 3815 : nfc->primaryVariable().dofIndicesNeighbor(),
2591 3815 : nfc->_connected_dof_indices,
2592 3815 : nfc->primaryVariable().scalingFactor());
2593 :
2594 : // We've handled Ken and Kne, finally handle Knn
2595 3815 : _fe_problem.cacheJacobianNeighbor(0);
2596 : }
2597 :
2598 : // Do the off-diagonals next
2599 3815 : const std::vector<MooseVariableFEBase *> coupled_vars = nfc->getCoupledMooseVars();
2600 7630 : for (const auto & jvar : coupled_vars)
2601 : {
2602 : // Only compute jacobians for nonlinear variables
2603 3815 : if (jvar->kind() != Moose::VAR_SOLVER)
2604 0 : continue;
2605 :
2606 : // Only compute Jacobian entries if this coupling is being used by the
2607 : // preconditioner
2608 3879 : if (nfc->variable().number() == jvar->number() ||
2609 128 : !_fe_problem.areCoupled(
2610 64 : nfc->variable().number(), jvar->number(), this->number()))
2611 3751 : continue;
2612 :
2613 : // Begin the off-diagonal node-face constraint accumulation phase for
2614 : // element and neighbor Jacobian blocks.
2615 64 : _fe_problem.prepareAssembly(0);
2616 64 : _fe_problem.prepareAssemblyNeighbor(0);
2617 :
2618 64 : nfc->prepareShapes(nfc->variable().number());
2619 64 : nfc->prepareNeighborShapes(jvar->number());
2620 :
2621 64 : nfc->computeOffDiagJacobian(jvar->number());
2622 :
2623 : // Cache the jacobian block for the secondary side
2624 128 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2625 64 : nfc->_Kee,
2626 : secondary_dofs,
2627 64 : nfc->_connected_dof_indices,
2628 : scaling_factor);
2629 :
2630 : // Cache Ken, Kne, Knn
2631 64 : if (nfc->addCouplingEntriesToJacobian())
2632 : {
2633 : // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2634 : // factor when we're overwriting secondary stuff)
2635 128 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2636 64 : nfc->_Ken,
2637 : secondary_dofs,
2638 64 : jvar->dofIndicesNeighbor(),
2639 : scaling_factor);
2640 :
2641 : // Use _connected_dof_indices to get all the correct columns
2642 128 : nfc->addJacobian(_fe_problem.assembly(0, number()),
2643 64 : nfc->_Kne,
2644 64 : nfc->variable().dofIndicesNeighbor(),
2645 64 : nfc->_connected_dof_indices,
2646 64 : nfc->variable().scalingFactor());
2647 :
2648 : // We've handled Ken and Kne, finally handle Knn
2649 64 : _fe_problem.cacheJacobianNeighbor(0);
2650 : }
2651 : }
2652 3815 : }
2653 : }
2654 : }
2655 : }
2656 : }
2657 : }
2658 1921 : if (_assemble_constraints_separately)
2659 : {
2660 : // See if constraints were applied anywhere
2661 0 : _communicator.max(constraints_applied);
2662 :
2663 0 : if (constraints_applied)
2664 : {
2665 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2666 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2667 : PETSC_TRUE));
2668 :
2669 0 : jacobian.close();
2670 0 : jacobian.zero_rows(zero_rows, 0.0);
2671 0 : jacobian.close();
2672 0 : _fe_problem.addCachedJacobian(0);
2673 0 : jacobian.close();
2674 : }
2675 : }
2676 : }
2677 4479 : if (!_assemble_constraints_separately)
2678 : {
2679 : // See if constraints were applied anywhere
2680 4479 : _communicator.max(constraints_applied);
2681 :
2682 4479 : if (constraints_applied)
2683 : {
2684 993 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2685 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2686 : PETSC_TRUE));
2687 :
2688 993 : jacobian.close();
2689 993 : jacobian.zero_rows(zero_rows, 0.0);
2690 993 : jacobian.close();
2691 993 : _fe_problem.addCachedJacobian(0);
2692 993 : jacobian.close();
2693 : }
2694 : }
2695 :
2696 4479 : THREAD_ID tid = 0;
2697 : // go over element-element constraint interface
2698 4479 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
2699 4479 : for (const auto & it : element_pair_locators)
2700 : {
2701 0 : ElementPairLocator & elem_pair_loc = *(it.second);
2702 :
2703 0 : if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
2704 : {
2705 : // ElemElemConstraint objects
2706 : const auto & element_constraints =
2707 0 : _constraints.getActiveElemElemConstraints(it.first, displaced);
2708 :
2709 : // go over pair elements
2710 : const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
2711 0 : elem_pair_loc.getElemPairs();
2712 0 : for (const auto & pr : elem_pairs)
2713 : {
2714 0 : const Elem * elem1 = pr.first;
2715 0 : const Elem * elem2 = pr.second;
2716 :
2717 0 : if (elem1->processor_id() != processor_id())
2718 0 : continue;
2719 :
2720 0 : const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
2721 :
2722 : // for each element process constraints on the
2723 0 : for (const auto & ec : element_constraints)
2724 : {
2725 0 : _fe_problem.setCurrentSubdomainID(elem1, tid);
2726 0 : subproblem.reinitElemPhys(elem1, info._elem1_constraint_q_point, tid);
2727 0 : _fe_problem.setNeighborSubdomainID(elem2, tid);
2728 0 : subproblem.reinitNeighborPhys(elem2, info._elem2_constraint_q_point, tid);
2729 :
2730 : // Begin the element-element constraint accumulation phase for element and neighbor
2731 : // Jacobian blocks.
2732 0 : _fe_problem.prepareAssembly(tid);
2733 0 : _fe_problem.prepareAssemblyNeighbor(tid);
2734 :
2735 0 : ec->prepareShapes(ec->variable().number());
2736 0 : ec->prepareNeighborShapes(ec->variable().number());
2737 :
2738 0 : ec->reinit(info);
2739 0 : ec->computeJacobian();
2740 0 : _fe_problem.cacheJacobian(tid);
2741 0 : _fe_problem.cacheJacobianNeighbor(tid);
2742 : }
2743 0 : _fe_problem.addCachedJacobian(tid);
2744 : }
2745 : }
2746 : }
2747 :
2748 : // go over NodeElemConstraints
2749 4479 : std::set<dof_id_type> unique_secondary_node_ids;
2750 4479 : constraints_applied = false;
2751 19139 : for (const auto & secondary_id : _mesh.meshSubdomains())
2752 : {
2753 72082 : for (const auto & primary_id : _mesh.meshSubdomains())
2754 : {
2755 57422 : if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
2756 : {
2757 : const auto & constraints =
2758 162 : _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
2759 :
2760 : // get unique set of ids of all nodes on current block
2761 162 : unique_secondary_node_ids.clear();
2762 162 : const MeshBase & meshhelper = _mesh.getMesh();
2763 324 : for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
2764 8046 : meshhelper.active_subdomain_elements_end(secondary_id)))
2765 : {
2766 44226 : for (auto & n : elem->node_ref_range())
2767 36666 : unique_secondary_node_ids.insert(n.id());
2768 162 : }
2769 :
2770 12636 : for (auto secondary_node_id : unique_secondary_node_ids)
2771 : {
2772 12474 : const Node & secondary_node = _mesh.nodeRef(secondary_node_id);
2773 : // check if secondary node is on current processor
2774 12474 : if (secondary_node.processor_id() == processor_id())
2775 : {
2776 : // This reinits the variables that exist on the secondary node
2777 9702 : _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
2778 :
2779 9702 : _fe_problem.reinitOffDiagScalars(0);
2780 :
2781 20636 : for (const auto & nec : constraints)
2782 : {
2783 10934 : if (nec->shouldApply())
2784 : {
2785 4970 : constraints_applied = true;
2786 :
2787 : // Begin the diagonal node-element constraint accumulation phase for
2788 : // element and neighbor Jacobian blocks.
2789 4970 : _fe_problem.prepareAssembly(0);
2790 4970 : _fe_problem.prepareAssemblyNeighbor(0);
2791 :
2792 4970 : nec->_jacobian = &jacobian_to_view;
2793 4970 : nec->prepareShapes(nec->variable().number());
2794 4970 : nec->prepareNeighborShapes(nec->variable().number());
2795 :
2796 4970 : nec->computeJacobian();
2797 :
2798 4970 : if (nec->overwriteSecondaryJacobian())
2799 : {
2800 : // Add this variable's dof's row to be zeroed
2801 0 : zero_rows.push_back(nec->variable().nodalDofIndex());
2802 : }
2803 :
2804 4970 : std::vector<dof_id_type> secondary_dofs(1, nec->variable().nodalDofIndex());
2805 :
2806 : // Cache the jacobian block for the secondary side
2807 9940 : nec->addJacobian(_fe_problem.assembly(0, number()),
2808 4970 : nec->_Kee,
2809 : secondary_dofs,
2810 4970 : nec->_connected_dof_indices,
2811 4970 : nec->variable().scalingFactor());
2812 :
2813 : // Cache the jacobian block for the primary side
2814 9940 : nec->addJacobian(_fe_problem.assembly(0, number()),
2815 4970 : nec->_Kne,
2816 4970 : nec->primaryVariable().dofIndicesNeighbor(),
2817 4970 : nec->_connected_dof_indices,
2818 4970 : nec->primaryVariable().scalingFactor());
2819 :
2820 4970 : _fe_problem.cacheJacobian(0);
2821 4970 : _fe_problem.cacheJacobianNeighbor(0);
2822 :
2823 : // Do the off-diagonals next
2824 4970 : const std::vector<MooseVariableFEBase *> coupled_vars = nec->getCoupledMooseVars();
2825 10010 : for (const auto & jvar : coupled_vars)
2826 : {
2827 : // Only compute jacobians for nonlinear variables
2828 5040 : if (jvar->kind() != Moose::VAR_SOLVER)
2829 70 : continue;
2830 :
2831 : // Only compute Jacobian entries if this coupling is being used by the
2832 : // preconditioner
2833 5530 : if (nec->variable().number() == jvar->number() ||
2834 1120 : !_fe_problem.areCoupled(
2835 560 : nec->variable().number(), jvar->number(), this->number()))
2836 4410 : continue;
2837 :
2838 : // Begin the off-diagonal node-element constraint accumulation phase for
2839 : // element and neighbor Jacobian blocks.
2840 560 : _fe_problem.prepareAssembly(0);
2841 560 : _fe_problem.prepareAssemblyNeighbor(0);
2842 :
2843 560 : nec->prepareShapes(nec->variable().number());
2844 560 : nec->prepareNeighborShapes(jvar->number());
2845 :
2846 560 : nec->computeOffDiagJacobian(jvar->number());
2847 :
2848 : // Cache the jacobian block for the secondary side
2849 1120 : nec->addJacobian(_fe_problem.assembly(0, number()),
2850 560 : nec->_Kee,
2851 : secondary_dofs,
2852 560 : nec->_connected_dof_indices,
2853 560 : nec->variable().scalingFactor());
2854 :
2855 : // Cache the jacobian block for the primary side
2856 1120 : nec->addJacobian(_fe_problem.assembly(0, number()),
2857 560 : nec->_Kne,
2858 560 : nec->variable().dofIndicesNeighbor(),
2859 560 : nec->_connected_dof_indices,
2860 560 : nec->variable().scalingFactor());
2861 :
2862 560 : _fe_problem.cacheJacobian(0);
2863 560 : _fe_problem.cacheJacobianNeighbor(0);
2864 : }
2865 4970 : }
2866 : }
2867 : }
2868 : }
2869 : }
2870 : }
2871 : }
2872 : // See if constraints were applied anywhere
2873 4479 : _communicator.max(constraints_applied);
2874 :
2875 4479 : if (constraints_applied)
2876 : {
2877 162 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2878 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2879 : PETSC_TRUE));
2880 :
2881 162 : jacobian.close();
2882 162 : jacobian.zero_rows(zero_rows, 0.0);
2883 162 : jacobian.close();
2884 162 : _fe_problem.addCachedJacobian(0);
2885 162 : jacobian.close();
2886 : }
2887 4479 : }
2888 :
2889 : void
2890 476278 : NonlinearSystemBase::computeScalarKernelsJacobians(const std::set<TagID> & tags)
2891 : {
2892 : MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
2893 :
2894 476278 : if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2895 467030 : scalar_kernel_warehouse = &_scalar_kernels;
2896 9248 : else if (tags.size() == 1)
2897 7813 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2898 : else
2899 1435 : scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagsObjectWarehouse(tags, 0));
2900 :
2901 : // Compute the diagonal block for scalar variables
2902 476278 : if (scalar_kernel_warehouse->hasActiveObjects())
2903 : {
2904 13504 : const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
2905 :
2906 13504 : _fe_problem.reinitScalars(/*tid=*/0);
2907 :
2908 13504 : _fe_problem.reinitOffDiagScalars(/*_tid*/ 0);
2909 :
2910 13504 : bool have_scalar_contributions = false;
2911 49561 : for (const auto & kernel : scalars)
2912 : {
2913 36057 : if (!kernel->computesJacobian())
2914 0 : continue;
2915 :
2916 36057 : kernel->reinit();
2917 36057 : const std::vector<dof_id_type> & dof_indices = kernel->variable().dofIndices();
2918 36057 : const DofMap & dof_map = kernel->variable().dofMap();
2919 36057 : const dof_id_type first_dof = dof_map.first_dof();
2920 36057 : const dof_id_type end_dof = dof_map.end_dof();
2921 42154 : for (dof_id_type dof : dof_indices)
2922 : {
2923 36099 : if (dof >= first_dof && dof < end_dof)
2924 : {
2925 30002 : kernel->computeJacobian();
2926 30002 : _fe_problem.addJacobianOffDiagScalar(kernel->variable().number());
2927 30002 : have_scalar_contributions = true;
2928 30002 : break;
2929 : }
2930 : }
2931 : }
2932 :
2933 13504 : if (have_scalar_contributions)
2934 11603 : _fe_problem.addJacobianScalar();
2935 : }
2936 476278 : }
2937 :
2938 : void
2939 491811 : NonlinearSystemBase::jacobianSetup()
2940 : {
2941 491811 : SolverSystem::jacobianSetup();
2942 :
2943 1035428 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2944 : {
2945 543617 : _kernels.jacobianSetup(tid);
2946 543617 : _nodal_kernels.jacobianSetup(tid);
2947 543617 : _dirac_kernels.jacobianSetup(tid);
2948 543617 : if (_doing_dg)
2949 2706 : _dg_kernels.jacobianSetup(tid);
2950 543617 : _interface_kernels.jacobianSetup(tid);
2951 543617 : _element_dampers.jacobianSetup(tid);
2952 543617 : _nodal_dampers.jacobianSetup(tid);
2953 543617 : _integrated_bcs.jacobianSetup(tid);
2954 : }
2955 491811 : _scalar_kernels.jacobianSetup();
2956 491811 : _constraints.jacobianSetup();
2957 491811 : _general_dampers.jacobianSetup();
2958 491811 : _nodal_bcs.jacobianSetup();
2959 491811 : _preset_nodal_bcs.jacobianSetup();
2960 491811 : _ad_preset_nodal_bcs.jacobianSetup();
2961 :
2962 : #ifdef MOOSE_KOKKOS_ENABLED
2963 358725 : _kokkos_kernels.jacobianSetup();
2964 358725 : _kokkos_nodal_kernels.jacobianSetup();
2965 358725 : _kokkos_integrated_bcs.jacobianSetup();
2966 358725 : _kokkos_nodal_bcs.jacobianSetup();
2967 : #endif
2968 :
2969 : // Avoid recursion
2970 491811 : if (this == &_fe_problem.currentNonlinearSystem())
2971 476278 : _fe_problem.jacobianSetup();
2972 491811 : }
2973 :
2974 : void
2975 476278 : NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
2976 : {
2977 1428834 : TIME_SECTION("computeJacobianInternal", 3);
2978 :
2979 476278 : _fe_problem.setCurrentNonlinearSystem(number());
2980 :
2981 : // Make matrix ready to use
2982 476278 : activateAllMatrixTags();
2983 :
2984 1421516 : for (auto tag : tags)
2985 : {
2986 945238 : if (!hasMatrix(tag))
2987 467030 : continue;
2988 :
2989 478208 : auto & jacobian = getMatrix(tag);
2990 : // Necessary for speed
2991 478208 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2992 : {
2993 477212 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2994 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2995 : PETSC_TRUE));
2996 477212 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2997 11991 : LibmeshPetscCall(
2998 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2999 477212 : if (_fe_problem.ignoreZerosInJacobian())
3000 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3001 : MAT_IGNORE_ZERO_ENTRIES,
3002 : PETSC_TRUE));
3003 : }
3004 : }
3005 :
3006 476278 : jacobianSetup();
3007 :
3008 : // Jacobian contributions from UOs - for now this is used for ray tracing
3009 : // and ray kernels that contribute to the Jacobian (think line sources)
3010 476278 : std::vector<UserObject *> uos;
3011 476278 : _fe_problem.theWarehouse()
3012 952556 : .query()
3013 476278 : .condition<AttribSystem>("UserObject")
3014 476278 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
3015 476278 : .queryInto(uos);
3016 476278 : for (auto & uo : uos)
3017 0 : uo->jacobianSetup();
3018 476278 : for (auto & uo : uos)
3019 : {
3020 0 : uo->initialize();
3021 0 : uo->execute();
3022 0 : uo->finalize();
3023 : }
3024 :
3025 : // reinit scalar variables
3026 1002916 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3027 526638 : _fe_problem.reinitScalars(tid);
3028 :
3029 : #ifdef MOOSE_KOKKOS_ENABLED
3030 347491 : if (_fe_problem.hasKokkosResidualObjects())
3031 10381 : computeKokkosJacobian(tags);
3032 : #endif
3033 :
3034 : PARALLEL_TRY
3035 : {
3036 : // We would like to compute ScalarKernels, block NodalKernels, FVFluxKernels, and mortar objects
3037 : // up front because we want these included whether we are computing an ordinary Jacobian or a
3038 : // Jacobian for determining variable scaling factors
3039 476278 : computeScalarKernelsJacobians(tags);
3040 :
3041 : // Block restricted Nodal Kernels
3042 476278 : if (_nodal_kernels.hasActiveBlockObjects())
3043 : {
3044 3698 : ComputeNodalKernelJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
3045 3698 : const ConstNodeRange & range = _fe_problem.getCurrentAlgebraicNodeRange();
3046 3698 : Threads::parallel_reduce(range, cnkjt);
3047 :
3048 3698 : unsigned int n_threads = libMesh::n_threads();
3049 8590 : for (unsigned int i = 0; i < n_threads;
3050 : i++) // Add any cached jacobians that might be hanging around
3051 4892 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
3052 3698 : }
3053 :
3054 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
3055 476278 : if (_fe_problem.haveFV())
3056 : {
3057 : // the same loop works for both residual and jacobians because it keys
3058 : // off of FEProblem's _currently_computing_jacobian parameter
3059 : ComputeFVFluxJacobianThread<FVRange> fvj(
3060 24156 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
3061 24156 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
3062 24156 : Threads::parallel_reduce(faces, fvj);
3063 24156 : }
3064 476278 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
3065 476278 : displaced_problem && displaced_problem->haveFV())
3066 : {
3067 : ComputeFVFluxJacobianThread<FVRange> fvr(
3068 0 : _fe_problem, this->number(), tags, /*on_displaced=*/true);
3069 0 : FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
3070 0 : displaced_problem->mesh().ownedFaceInfoEnd());
3071 0 : Threads::parallel_reduce(faces, fvr);
3072 476278 : }
3073 :
3074 476278 : mortarConstraints(Moose::ComputeType::Jacobian, {}, tags);
3075 :
3076 : // Get our element range for looping over
3077 476278 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
3078 :
3079 476278 : if (_fe_problem.computingScalingJacobian())
3080 : {
3081 : // Only compute Jacobians corresponding to the diagonals of volumetric compute objects
3082 : // because this typically gives us a good representation of the physics. NodalBCs and
3083 : // Constraints can introduce dramatically different scales (often order unity).
3084 : // IntegratedBCs and/or InterfaceKernels may use penalty factors. DGKernels may be ok, but
3085 : // they are almost always used in conjunction with Kernels
3086 565 : ComputeJacobianForScalingThread cj(_fe_problem, tags);
3087 565 : Threads::parallel_reduce(elem_range, cj);
3088 565 : unsigned int n_threads = libMesh::n_threads();
3089 1200 : for (unsigned int i = 0; i < n_threads;
3090 : i++) // Add any Jacobian contributions still hanging around
3091 635 : _fe_problem.addCachedJacobian(i);
3092 :
3093 : // Check whether any exceptions were thrown and propagate this information for parallel
3094 : // consistency before
3095 : // 1) we do parallel communication when closing tagged matrices
3096 : // 2) early returning before reaching our PARALLEL_CATCH below
3097 565 : _fe_problem.checkExceptionAndStopSolve();
3098 :
3099 565 : closeTaggedMatrices(tags);
3100 :
3101 565 : return;
3102 565 : }
3103 :
3104 475713 : switch (_fe_problem.coupling())
3105 : {
3106 397098 : case Moose::COUPLING_DIAG:
3107 : {
3108 397098 : ComputeJacobianThread cj(_fe_problem, tags);
3109 397098 : Threads::parallel_reduce(elem_range, cj);
3110 :
3111 397094 : unsigned int n_threads = libMesh::n_threads();
3112 837008 : for (unsigned int i = 0; i < n_threads;
3113 : i++) // Add any Jacobian contributions still hanging around
3114 439914 : _fe_problem.addCachedJacobian(i);
3115 :
3116 : // Boundary restricted Nodal Kernels
3117 397094 : if (_nodal_kernels.hasActiveBoundaryObjects())
3118 : {
3119 42 : ComputeNodalKernelBCJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
3120 42 : const ConstBndNodeRange & bnd_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
3121 :
3122 42 : Threads::parallel_reduce(bnd_range, cnkjt);
3123 42 : unsigned int n_threads = libMesh::n_threads();
3124 88 : for (unsigned int i = 0; i < n_threads;
3125 : i++) // Add any cached jacobians that might be hanging around
3126 46 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
3127 42 : }
3128 397094 : }
3129 397094 : break;
3130 :
3131 78615 : default:
3132 : case Moose::COUPLING_CUSTOM:
3133 : {
3134 78615 : ComputeFullJacobianThread cj(_fe_problem, tags);
3135 78615 : Threads::parallel_reduce(elem_range, cj);
3136 78615 : unsigned int n_threads = libMesh::n_threads();
3137 :
3138 164696 : for (unsigned int i = 0; i < n_threads; i++)
3139 86084 : _fe_problem.addCachedJacobian(i);
3140 :
3141 : // Boundary restricted Nodal Kernels
3142 78612 : if (_nodal_kernels.hasActiveBoundaryObjects())
3143 : {
3144 9 : ComputeNodalKernelBCJacobiansThread cnkjt(_fe_problem, *this, _nodal_kernels, tags);
3145 9 : const ConstBndNodeRange & bnd_range = _fe_problem.getCurrentAlgebraicBndNodeRange();
3146 :
3147 9 : Threads::parallel_reduce(bnd_range, cnkjt);
3148 9 : unsigned int n_threads = libMesh::n_threads();
3149 19 : for (unsigned int i = 0; i < n_threads;
3150 : i++) // Add any cached jacobians that might be hanging around
3151 10 : _fe_problem.assembly(i, number()).addCachedJacobian(Assembly::GlobalDataKey{});
3152 9 : }
3153 78615 : }
3154 78612 : break;
3155 : }
3156 :
3157 475706 : computeDiracContributions(tags, true);
3158 :
3159 : static bool first = true;
3160 :
3161 : // This adds zeroes into geometric coupling entries to ensure they stay in the matrix
3162 950834 : if ((_fe_problem.restoreOriginalNonzeroPattern() || first) &&
3163 475131 : _add_implicit_geometric_coupling_entries_to_jacobian)
3164 : {
3165 769 : first = false;
3166 769 : addImplicitGeometricCouplingEntries(_fe_problem.geomSearchData());
3167 :
3168 769 : if (_fe_problem.getDisplacedProblem())
3169 18 : addImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData());
3170 : }
3171 : }
3172 475703 : PARALLEL_CATCH;
3173 :
3174 : // Have no idea how to have constraints work
3175 : // with the tag system
3176 : PARALLEL_TRY
3177 : {
3178 : // Add in Jacobian contributions from other Constraints
3179 475640 : if (_fe_problem._has_constraints && tags.count(systemMatrixTag()))
3180 : {
3181 : // Some constraints need to be able to read values from the Jacobian, which requires that it
3182 : // be closed/assembled
3183 3476 : auto & system_matrix = getMatrix(systemMatrixTag());
3184 3476 : std::unique_ptr<SparseMatrix<Number>> hash_copy;
3185 : const SparseMatrix<Number> * view_jac_ptr;
3186 3627 : auto make_readable_jacobian = [&]()
3187 : {
3188 : #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
3189 3627 : if (system_matrix.use_hash_table())
3190 : {
3191 2169 : hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
3192 2169 : view_jac_ptr = hash_copy.get();
3193 : }
3194 : else
3195 1458 : view_jac_ptr = &system_matrix;
3196 : #else
3197 : view_jac_ptr = &system_matrix;
3198 : #endif
3199 3627 : if (view_jac_ptr == &system_matrix)
3200 1458 : system_matrix.close();
3201 3627 : };
3202 :
3203 3476 : make_readable_jacobian();
3204 :
3205 : // Nodal Constraints
3206 3476 : const bool had_nodal_constraints = enforceNodalConstraintsJacobian(*view_jac_ptr);
3207 3476 : if (had_nodal_constraints)
3208 : // We have to make a new readable Jacobian
3209 151 : make_readable_jacobian();
3210 :
3211 : // Undisplaced Constraints
3212 3476 : constraintJacobians(*view_jac_ptr, false);
3213 :
3214 : // Displaced Constraints
3215 3476 : if (_fe_problem.getDisplacedProblem())
3216 1003 : constraintJacobians(*view_jac_ptr, true);
3217 3476 : }
3218 : }
3219 475640 : PARALLEL_CATCH;
3220 :
3221 475640 : computeNodalBCsJacobian(tags);
3222 475640 : closeTaggedMatrices(tags);
3223 :
3224 : // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
3225 : // on boundary nodes
3226 475640 : if (_has_nodalbc_diag_save_in)
3227 21 : _fe_problem.getAuxiliarySystem().solution().close();
3228 :
3229 475640 : if (hasDiagSaveIn())
3230 170 : _fe_problem.getAuxiliarySystem().update();
3231 :
3232 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
3233 : // counters
3234 475640 : _app.solutionInvalidity().syncIteration();
3235 475640 : _app.solutionInvalidity().accumulateIterationIntoTimeStepOccurences();
3236 476902 : }
3237 :
3238 : void
3239 0 : NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian)
3240 : {
3241 0 : _nl_matrix_tags.clear();
3242 :
3243 0 : auto & tags = _fe_problem.getMatrixTags();
3244 :
3245 0 : for (auto & tag : tags)
3246 0 : _nl_matrix_tags.insert(tag.second);
3247 :
3248 0 : computeJacobian(jacobian, _nl_matrix_tags);
3249 0 : }
3250 :
3251 : void
3252 0 : NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian, const std::set<TagID> & tags)
3253 : {
3254 0 : associateMatrixToTag(jacobian, systemMatrixTag());
3255 :
3256 0 : computeJacobianTags(tags);
3257 :
3258 0 : disassociateMatrixFromTag(jacobian, systemMatrixTag());
3259 0 : }
3260 :
3261 : void
3262 476278 : NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
3263 : {
3264 1428834 : TIME_SECTION("computeJacobianTags", 5);
3265 :
3266 476278 : FloatingPointExceptionGuard fpe_guard(_app);
3267 :
3268 : try
3269 : {
3270 476278 : computeJacobianInternal(tags);
3271 : }
3272 66 : catch (MooseException & e)
3273 : {
3274 : // The buck stops here, we have already handled the exception by
3275 : // calling stopSolve(), it is now up to PETSc to return a
3276 : // "diverged" reason during the next solve.
3277 63 : }
3278 476274 : }
3279 :
3280 : void
3281 263 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
3282 : {
3283 263 : _nl_matrix_tags.clear();
3284 :
3285 263 : auto & tags = _fe_problem.getMatrixTags();
3286 789 : for (auto & tag : tags)
3287 526 : _nl_matrix_tags.insert(tag.second);
3288 :
3289 263 : computeJacobianBlocks(blocks, _nl_matrix_tags);
3290 263 : }
3291 :
3292 : void
3293 587 : NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
3294 : const std::set<TagID> & tags)
3295 : {
3296 1761 : TIME_SECTION("computeJacobianBlocks", 3);
3297 587 : FloatingPointExceptionGuard fpe_guard(_app);
3298 :
3299 2050 : for (unsigned int i = 0; i < blocks.size(); i++)
3300 : {
3301 1463 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3302 :
3303 1463 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3304 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
3305 : PETSC_TRUE));
3306 1463 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
3307 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3308 : MAT_NEW_NONZERO_ALLOCATION_ERR,
3309 : PETSC_TRUE));
3310 :
3311 1463 : jacobian.zero();
3312 : }
3313 :
3314 1237 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3315 650 : _fe_problem.reinitScalars(tid);
3316 :
3317 : PARALLEL_TRY
3318 : {
3319 587 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
3320 587 : ComputeJacobianBlocksThread cjb(_fe_problem, blocks, tags);
3321 587 : Threads::parallel_reduce(elem_range, cjb);
3322 587 : }
3323 587 : PARALLEL_CATCH;
3324 :
3325 2050 : for (unsigned int i = 0; i < blocks.size(); i++)
3326 1463 : blocks[i]->_jacobian.close();
3327 :
3328 2050 : for (unsigned int i = 0; i < blocks.size(); i++)
3329 : {
3330 1463 : libMesh::System & precond_system = blocks[i]->_precond_system;
3331 1463 : SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3332 :
3333 1463 : unsigned int ivar = blocks[i]->_ivar;
3334 1463 : unsigned int jvar = blocks[i]->_jvar;
3335 :
3336 : // Dirichlet BCs
3337 1463 : std::vector<numeric_index_type> zero_rows;
3338 : PARALLEL_TRY
3339 : {
3340 1463 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3341 49739 : for (const auto & bnode : bnd_nodes)
3342 : {
3343 48276 : BoundaryID boundary_id = bnode->_bnd_id;
3344 48276 : Node * node = bnode->_node;
3345 :
3346 48276 : if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
3347 : {
3348 42756 : const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
3349 :
3350 42756 : if (node->processor_id() == processor_id())
3351 : {
3352 32552 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3353 :
3354 150180 : for (const auto & bc : bcs)
3355 117628 : if (bc->variable().number() == ivar && bc->shouldApply())
3356 : {
3357 : // The first zero is for the variable number... there is only one variable in
3358 : // each mini-system The second zero only works with Lagrange elements!
3359 44344 : zero_rows.push_back(node->dof_number(precond_system.number(), 0, 0));
3360 : }
3361 : }
3362 : }
3363 : }
3364 : }
3365 1463 : PARALLEL_CATCH;
3366 :
3367 1463 : jacobian.close();
3368 :
3369 : // This zeroes the rows corresponding to Dirichlet BCs and puts a 1.0 on the diagonal
3370 1463 : if (ivar == jvar)
3371 1366 : jacobian.zero_rows(zero_rows, 1.0);
3372 : else
3373 97 : jacobian.zero_rows(zero_rows, 0.0);
3374 :
3375 1463 : jacobian.close();
3376 1463 : }
3377 587 : }
3378 :
3379 : void
3380 346870 : NonlinearSystemBase::updateActive(THREAD_ID tid)
3381 : {
3382 346870 : _element_dampers.updateActive(tid);
3383 346870 : _nodal_dampers.updateActive(tid);
3384 346870 : _integrated_bcs.updateActive(tid);
3385 346870 : _dg_kernels.updateActive(tid);
3386 346870 : _interface_kernels.updateActive(tid);
3387 346870 : _dirac_kernels.updateActive(tid);
3388 346870 : _kernels.updateActive(tid);
3389 346870 : _nodal_kernels.updateActive(tid);
3390 :
3391 346870 : if (tid == 0)
3392 : {
3393 315532 : _general_dampers.updateActive();
3394 315532 : _nodal_bcs.updateActive();
3395 315532 : _preset_nodal_bcs.updateActive();
3396 315532 : _ad_preset_nodal_bcs.updateActive();
3397 315532 : _splits.updateActive();
3398 315532 : _constraints.updateActive();
3399 315532 : _scalar_kernels.updateActive();
3400 :
3401 : #ifdef MOOSE_KOKKOS_ENABLED
3402 230284 : _kokkos_kernels.updateActive();
3403 230284 : _kokkos_nodal_kernels.updateActive();
3404 230284 : _kokkos_integrated_bcs.updateActive();
3405 230284 : _kokkos_nodal_bcs.updateActive();
3406 230284 : _kokkos_preset_nodal_bcs.updateActive();
3407 : #endif
3408 : }
3409 346870 : }
3410 :
3411 : Real
3412 1625 : NonlinearSystemBase::computeDamping(const NumericVector<Number> & solution,
3413 : const NumericVector<Number> & update)
3414 : {
3415 : // Default to no damping
3416 1625 : Real damping = 1.0;
3417 1625 : bool has_active_dampers = false;
3418 :
3419 : try
3420 : {
3421 1625 : if (_element_dampers.hasActiveObjects())
3422 : {
3423 : PARALLEL_TRY
3424 : {
3425 2400 : TIME_SECTION("computeDampers", 3, "Computing Dampers");
3426 480 : has_active_dampers = true;
3427 480 : *_increment_vec = update;
3428 480 : ComputeElemDampingThread cid(_fe_problem, *this);
3429 480 : Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicElementRange(), cid);
3430 480 : damping = std::min(cid.damping(), damping);
3431 480 : }
3432 480 : PARALLEL_CATCH;
3433 : }
3434 :
3435 1560 : if (_nodal_dampers.hasActiveObjects())
3436 : {
3437 : PARALLEL_TRY
3438 : {
3439 3855 : TIME_SECTION("computeDamping::element", 3, "Computing Element Damping");
3440 :
3441 771 : has_active_dampers = true;
3442 771 : *_increment_vec = update;
3443 771 : ComputeNodalDampingThread cndt(_fe_problem, *this);
3444 771 : Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicNodeRange(), cndt);
3445 771 : damping = std::min(cndt.damping(), damping);
3446 771 : }
3447 771 : PARALLEL_CATCH;
3448 : }
3449 :
3450 1501 : if (_general_dampers.hasActiveObjects())
3451 : {
3452 : PARALLEL_TRY
3453 : {
3454 1830 : TIME_SECTION("computeDamping::general", 3, "Computing General Damping");
3455 :
3456 366 : has_active_dampers = true;
3457 366 : const auto & gdampers = _general_dampers.getActiveObjects();
3458 732 : for (const auto & damper : gdampers)
3459 : {
3460 366 : Real gd_damping = damper->computeDamping(solution, update);
3461 : try
3462 : {
3463 366 : damper->checkMinDamping(gd_damping);
3464 : }
3465 6 : catch (MooseException & e)
3466 : {
3467 12 : _fe_problem.setException(e.what());
3468 6 : }
3469 366 : damping = std::min(gd_damping, damping);
3470 : }
3471 366 : }
3472 366 : PARALLEL_CATCH;
3473 : }
3474 : }
3475 130 : catch (MooseException & e)
3476 : {
3477 : // The buck stops here, we have already handled the exception by
3478 : // calling stopSolve(), it is now up to PETSc to return a
3479 : // "diverged" reason during the next solve.
3480 130 : }
3481 0 : catch (std::exception & e)
3482 : {
3483 : // Allow the libmesh error/exception on negative jacobian
3484 0 : const std::string & message = e.what();
3485 0 : if (message.find("Jacobian") == std::string::npos)
3486 0 : throw;
3487 0 : }
3488 :
3489 1625 : _communicator.min(damping);
3490 :
3491 1625 : if (has_active_dampers && damping < 1.0)
3492 1221 : _console << " Damping factor: " << damping << std::endl;
3493 :
3494 1625 : return damping;
3495 : }
3496 :
3497 : void
3498 3530960 : NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
3499 : {
3500 3530960 : _fe_problem.clearDiracInfo();
3501 :
3502 3530960 : std::set<const Elem *> dirac_elements;
3503 :
3504 3530960 : if (_dirac_kernels.hasActiveObjects())
3505 : {
3506 177635 : TIME_SECTION("computeDirac", 3, "Computing DiracKernels");
3507 :
3508 : // TODO: Need a threading fix... but it's complicated!
3509 74159 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3510 : {
3511 38641 : const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
3512 89001 : for (const auto & dkernel : dkernels)
3513 : {
3514 50369 : dkernel->clearPoints();
3515 50369 : dkernel->addPoints();
3516 : }
3517 : }
3518 :
3519 35518 : ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
3520 :
3521 35518 : _fe_problem.getDiracElements(dirac_elements);
3522 :
3523 35518 : DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
3524 : // TODO: Make Dirac work thread!
3525 : // Threads::parallel_reduce(range, cd);
3526 :
3527 35518 : cd(range);
3528 :
3529 35518 : if (is_jacobian)
3530 5587 : for (const auto tid : make_range(libMesh::n_threads()))
3531 2921 : _fe_problem.addCachedJacobian(tid);
3532 35518 : }
3533 3530951 : }
3534 :
3535 : NumericVector<Number> &
3536 0 : NonlinearSystemBase::residualCopy()
3537 : {
3538 0 : if (!_residual_copy.get())
3539 0 : _residual_copy = NumericVector<Number>::build(_communicator);
3540 :
3541 0 : return *_residual_copy;
3542 : }
3543 :
3544 : NumericVector<Number> &
3545 344 : NonlinearSystemBase::residualGhosted()
3546 : {
3547 344 : _need_residual_ghosted = true;
3548 344 : if (!_residual_ghosted)
3549 : {
3550 : // The first time we realize we need a ghosted residual vector,
3551 : // we add it.
3552 556 : _residual_ghosted = &addVector("residual_ghosted", false, GHOSTED);
3553 :
3554 : // If we've already realized we need time and/or non-time
3555 : // residual vectors, but we haven't yet realized they need to be
3556 : // ghosted, fix that now.
3557 : //
3558 : // If an application changes its mind, the libMesh API lets us
3559 : // change the vector.
3560 278 : if (_Re_time)
3561 : {
3562 57 : const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
3563 57 : _Re_time = &system().add_vector(vector_name, false, GHOSTED);
3564 57 : }
3565 278 : if (_Re_non_time)
3566 : {
3567 278 : const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
3568 278 : _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
3569 278 : }
3570 : }
3571 344 : return *_residual_ghosted;
3572 : }
3573 :
3574 : void
3575 68567 : NonlinearSystemBase::augmentSparsity(SparsityPattern::Graph & sparsity,
3576 : std::vector<dof_id_type> & n_nz,
3577 : std::vector<dof_id_type> & n_oz)
3578 : {
3579 68567 : if (_add_implicit_geometric_coupling_entries_to_jacobian)
3580 : {
3581 1376 : _fe_problem.updateGeomSearch();
3582 :
3583 1376 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
3584 :
3585 1376 : findImplicitGeometricCouplingEntries(_fe_problem.geomSearchData(), graph);
3586 :
3587 1376 : if (_fe_problem.getDisplacedProblem())
3588 183 : findImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData(),
3589 : graph);
3590 :
3591 1376 : const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
3592 1376 : const dof_id_type end_dof_on_proc = dofMap().end_dof(processor_id());
3593 :
3594 : // The total number of dofs on and off processor
3595 1376 : const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
3596 1376 : const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
3597 :
3598 2338 : for (const auto & git : graph)
3599 : {
3600 962 : dof_id_type dof = git.first;
3601 962 : dof_id_type local_dof = dof - first_dof_on_proc;
3602 :
3603 962 : if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
3604 176 : continue;
3605 :
3606 786 : const auto & row = git.second;
3607 :
3608 786 : SparsityPattern::Row & sparsity_row = sparsity[local_dof];
3609 :
3610 786 : unsigned int original_row_length = sparsity_row.size();
3611 :
3612 786 : sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
3613 :
3614 1572 : SparsityPattern::sort_row(
3615 786 : sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
3616 :
3617 : // Fix up nonzero arrays
3618 5036 : for (const auto & coupled_dof : row)
3619 : {
3620 4250 : if (coupled_dof < first_dof_on_proc || coupled_dof >= end_dof_on_proc)
3621 : {
3622 1296 : if (n_oz[local_dof] < n_dofs_not_on_proc)
3623 648 : n_oz[local_dof]++;
3624 : }
3625 : else
3626 : {
3627 3602 : if (n_nz[local_dof] < n_dofs_on_proc)
3628 3602 : n_nz[local_dof]++;
3629 : }
3630 : }
3631 : }
3632 1376 : }
3633 68567 : }
3634 :
3635 : void
3636 0 : NonlinearSystemBase::setSolutionUDot(const NumericVector<Number> & u_dot)
3637 : {
3638 0 : *_u_dot = u_dot;
3639 0 : }
3640 :
3641 : void
3642 0 : NonlinearSystemBase::setSolutionUDotDot(const NumericVector<Number> & u_dotdot)
3643 : {
3644 0 : *_u_dotdot = u_dotdot;
3645 0 : }
3646 :
3647 : void
3648 0 : NonlinearSystemBase::setSolutionUDotOld(const NumericVector<Number> & u_dot_old)
3649 : {
3650 0 : *_u_dot_old = u_dot_old;
3651 0 : }
3652 :
3653 : void
3654 0 : NonlinearSystemBase::setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old)
3655 : {
3656 0 : *_u_dotdot_old = u_dotdot_old;
3657 0 : }
3658 :
3659 : void
3660 13557 : NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
3661 : {
3662 13557 : if (_preconditioner.get() != nullptr)
3663 4 : mooseError("More than one active Preconditioner detected");
3664 :
3665 13553 : _preconditioner = pc;
3666 13553 : }
3667 :
3668 : MoosePreconditioner const *
3669 55298 : NonlinearSystemBase::getPreconditioner() const
3670 : {
3671 55298 : return _preconditioner.get();
3672 : }
3673 :
3674 : void
3675 321 : NonlinearSystemBase::setupDampers()
3676 : {
3677 321 : _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
3678 321 : }
3679 :
3680 : void
3681 1480 : NonlinearSystemBase::reinitIncrementAtQpsForDampers(THREAD_ID /*tid*/,
3682 : const std::set<MooseVariable *> & damped_vars)
3683 : {
3684 2979 : for (const auto & var : damped_vars)
3685 1499 : var->computeIncrementAtQps(*_increment_vec);
3686 1480 : }
3687 :
3688 : void
3689 7473 : NonlinearSystemBase::reinitIncrementAtNodeForDampers(THREAD_ID /*tid*/,
3690 : const std::set<MooseVariable *> & damped_vars)
3691 : {
3692 14946 : for (const auto & var : damped_vars)
3693 7473 : var->computeIncrementAtNode(*_increment_vec);
3694 7473 : }
3695 :
3696 : void
3697 41035 : NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
3698 : {
3699 : // Obtain all blocks and variables covered by all kernels
3700 41035 : std::set<SubdomainID> input_subdomains;
3701 41035 : std::set<std::string> kernel_variables;
3702 :
3703 41035 : bool global_kernels_exist = false;
3704 41035 : global_kernels_exist |= _scalar_kernels.hasActiveObjects();
3705 41035 : global_kernels_exist |= _nodal_kernels.hasActiveObjects();
3706 :
3707 41035 : _kernels.subdomainsCovered(input_subdomains, kernel_variables);
3708 41035 : _dg_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3709 41035 : _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3710 41035 : _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3711 41035 : _constraints.subdomainsCovered(input_subdomains, kernel_variables);
3712 :
3713 : #ifdef MOOSE_KOKKOS_ENABLED
3714 30851 : _kokkos_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3715 30851 : _kokkos_nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3716 : #endif
3717 :
3718 41035 : if (_fe_problem.haveFV())
3719 : {
3720 2675 : std::vector<FVElementalKernel *> fv_elemental_kernels;
3721 2675 : _fe_problem.theWarehouse()
3722 5350 : .query()
3723 2675 : .template condition<AttribSystem>("FVElementalKernel")
3724 2675 : .queryInto(fv_elemental_kernels);
3725 :
3726 5596 : for (auto fv_kernel : fv_elemental_kernels)
3727 : {
3728 2921 : if (fv_kernel->blockRestricted())
3729 1334 : for (auto block_id : fv_kernel->blockIDs())
3730 733 : input_subdomains.insert(block_id);
3731 : else
3732 2320 : global_kernels_exist = true;
3733 2921 : kernel_variables.insert(fv_kernel->variable().name());
3734 :
3735 : // Check for lagrange multiplier
3736 2921 : if (dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel))
3737 248 : kernel_variables.insert(dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel)
3738 248 : ->lambdaVariable()
3739 124 : .name());
3740 : }
3741 :
3742 2675 : std::vector<FVFluxKernel *> fv_flux_kernels;
3743 2675 : _fe_problem.theWarehouse()
3744 5350 : .query()
3745 2675 : .template condition<AttribSystem>("FVFluxKernel")
3746 2675 : .queryInto(fv_flux_kernels);
3747 :
3748 6656 : for (auto fv_kernel : fv_flux_kernels)
3749 : {
3750 3981 : if (fv_kernel->blockRestricted())
3751 1688 : for (auto block_id : fv_kernel->blockIDs())
3752 892 : input_subdomains.insert(block_id);
3753 : else
3754 3185 : global_kernels_exist = true;
3755 3981 : kernel_variables.insert(fv_kernel->variable().name());
3756 : }
3757 :
3758 2675 : std::vector<FVInterfaceKernel *> fv_interface_kernels;
3759 2675 : _fe_problem.theWarehouse()
3760 5350 : .query()
3761 2675 : .template condition<AttribSystem>("FVInterfaceKernel")
3762 2675 : .queryInto(fv_interface_kernels);
3763 :
3764 2904 : for (auto fvik : fv_interface_kernels)
3765 229 : if (auto scalar_fvik = dynamic_cast<FVScalarLagrangeMultiplierInterface *>(fvik))
3766 13 : kernel_variables.insert(scalar_fvik->lambdaVariable().name());
3767 :
3768 2675 : std::vector<FVFluxBC *> fv_flux_bcs;
3769 2675 : _fe_problem.theWarehouse()
3770 5350 : .query()
3771 2675 : .template condition<AttribSystem>("FVFluxBC")
3772 2675 : .queryInto(fv_flux_bcs);
3773 :
3774 4026 : for (auto fvbc : fv_flux_bcs)
3775 1351 : if (auto scalar_fvbc = dynamic_cast<FVBoundaryScalarLagrangeMultiplierConstraint *>(fvbc))
3776 13 : kernel_variables.insert(scalar_fvbc->lambdaVariable().name());
3777 2675 : }
3778 :
3779 48274 : for (const auto & ibc : _integrated_bcs.getActiveObjects())
3780 : {
3781 7239 : const auto additional_variables_covered = ibc->additionalROVariables();
3782 7239 : kernel_variables.insert(additional_variables_covered.begin(),
3783 : additional_variables_covered.end());
3784 7239 : }
3785 :
3786 : // Check kernel coverage of subdomains (blocks) in your mesh
3787 41035 : if (!global_kernels_exist)
3788 : {
3789 37857 : std::set<SubdomainID> difference;
3790 37857 : std::set_difference(mesh_subdomains.begin(),
3791 : mesh_subdomains.end(),
3792 : input_subdomains.begin(),
3793 : input_subdomains.end(),
3794 : std::inserter(difference, difference.end()));
3795 :
3796 : // there supposed to be no kernels on this lower-dimensional subdomain
3797 38076 : for (const auto & id : _mesh.interiorLowerDBlocks())
3798 219 : difference.erase(id);
3799 38076 : for (const auto & id : _mesh.boundaryLowerDBlocks())
3800 219 : difference.erase(id);
3801 :
3802 37857 : if (!difference.empty())
3803 : {
3804 : std::vector<SubdomainID> difference_vec =
3805 6 : std::vector<SubdomainID>(difference.begin(), difference.end());
3806 6 : std::vector<SubdomainName> difference_names = _mesh.getSubdomainNames(difference_vec);
3807 6 : std::stringstream missing_block_names;
3808 6 : std::copy(difference_names.begin(),
3809 : difference_names.end(),
3810 6 : std::ostream_iterator<std::string>(missing_block_names, " "));
3811 6 : std::stringstream missing_block_ids;
3812 6 : std::copy(difference.begin(),
3813 : difference.end(),
3814 6 : std::ostream_iterator<unsigned int>(missing_block_ids, " "));
3815 :
3816 6 : mooseError("Each subdomain must contain at least one Kernel.\nThe following block(s) lack an "
3817 6 : "active kernel: " +
3818 6 : missing_block_names.str(),
3819 : " (ids: ",
3820 6 : missing_block_ids.str(),
3821 : ")");
3822 0 : }
3823 37851 : }
3824 :
3825 : // Check kernel use of variables
3826 41029 : std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
3827 :
3828 41029 : std::set<VariableName> difference;
3829 41029 : std::set_difference(variables.begin(),
3830 : variables.end(),
3831 : kernel_variables.begin(),
3832 : kernel_variables.end(),
3833 : std::inserter(difference, difference.end()));
3834 :
3835 : // skip checks for varaibles defined on lower-dimensional subdomain
3836 41029 : std::set<VariableName> vars(difference);
3837 41362 : for (auto & var_name : vars)
3838 : {
3839 333 : auto blks = getSubdomainsForVar(var_name);
3840 684 : for (const auto & id : blks)
3841 351 : if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
3842 351 : difference.erase(var_name);
3843 333 : }
3844 :
3845 41029 : if (!difference.empty())
3846 : {
3847 6 : std::stringstream missing_kernel_vars;
3848 6 : std::copy(difference.begin(),
3849 : difference.end(),
3850 6 : std::ostream_iterator<std::string>(missing_kernel_vars, " "));
3851 6 : mooseError("Each variable must be referenced by at least one active Kernel.\nThe following "
3852 6 : "variable(s) lack an active kernel: " +
3853 6 : missing_kernel_vars.str());
3854 0 : }
3855 41023 : }
3856 :
3857 : bool
3858 28448 : NonlinearSystemBase::containsTimeKernel()
3859 : {
3860 28448 : auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3861 :
3862 28448 : return time_kernels.hasActiveObjects();
3863 : }
3864 :
3865 : std::vector<std::string>
3866 0 : NonlinearSystemBase::timeKernelVariableNames()
3867 : {
3868 0 : std::vector<std::string> variable_names;
3869 0 : const auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3870 0 : if (time_kernels.hasActiveObjects())
3871 0 : for (const auto & kernel : time_kernels.getObjects())
3872 0 : variable_names.push_back(kernel->variable().name());
3873 :
3874 0 : return variable_names;
3875 0 : }
3876 :
3877 : bool
3878 27690 : NonlinearSystemBase::needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3879 : {
3880 : // IntegratedBCs are for now the only objects we consider to be consuming
3881 : // matprops on boundaries.
3882 27690 : if (_integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid))
3883 5409 : for (const auto & bc : _integrated_bcs.getActiveBoundaryObjects(bnd_id, tid))
3884 3278 : if (std::static_pointer_cast<MaterialPropertyInterface>(bc)->getMaterialPropertyCalled())
3885 974 : return true;
3886 :
3887 : // Thin layer heat transfer in the heat_transfer module is being used on a boundary even though
3888 : // it's an interface kernel. That boundary is external, on both sides of a gap in a mesh
3889 26716 : if (_interface_kernels.hasActiveBoundaryObjects(bnd_id, tid))
3890 321 : for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3891 252 : if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3892 183 : return true;
3893 :
3894 : // Because MortarConstraints do not inherit from BoundaryRestrictable, they are not sorted
3895 : // by boundary in the MooseObjectWarehouse. So for now, we return true for all boundaries
3896 : // Note: constraints are not threaded at this time
3897 26533 : if (_constraints.hasActiveObjects(/*tid*/ 0))
3898 2329 : for (const auto & ct : _constraints.getActiveObjects(/*tid*/ 0))
3899 1964 : if (auto mpi = std::dynamic_pointer_cast<MaterialPropertyInterface>(ct);
3900 1964 : mpi && mpi->getMaterialPropertyCalled())
3901 1964 : return true;
3902 25606 : return false;
3903 : }
3904 :
3905 : bool
3906 2700 : NonlinearSystemBase::needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3907 : {
3908 : // InterfaceKernels are for now the only objects we consider to be consuming matprops on internal
3909 : // boundaries.
3910 2700 : if (_interface_kernels.hasActiveBoundaryObjects(bnd_id, tid))
3911 398 : for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3912 307 : if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3913 216 : return true;
3914 2484 : return false;
3915 : }
3916 :
3917 : bool
3918 12223 : NonlinearSystemBase::needInternalNeighborSideMaterial(SubdomainID subdomain_id, THREAD_ID tid) const
3919 : {
3920 : // DGKernels are for now the only objects we consider to be consuming matprops on
3921 : // internal sides.
3922 12223 : if (_dg_kernels.hasActiveBlockObjects(subdomain_id, tid))
3923 677 : for (const auto & dg : _dg_kernels.getActiveBlockObjects(subdomain_id, tid))
3924 554 : if (std::static_pointer_cast<MaterialPropertyInterface>(dg)->getMaterialPropertyCalled())
3925 410 : return true;
3926 : // NOTE:
3927 : // HDG kernels do not require face material properties on internal sides at this time.
3928 : // The idea is to have element locality of HDG for hybridization
3929 11813 : return false;
3930 : }
3931 :
3932 : bool
3933 0 : NonlinearSystemBase::doingDG() const
3934 : {
3935 0 : return _doing_dg;
3936 : }
3937 :
3938 : void
3939 503 : NonlinearSystemBase::setPreviousNewtonSolution(const NumericVector<Number> & soln)
3940 : {
3941 503 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
3942 503 : getVector(Moose::PREVIOUS_NL_SOLUTION_TAG) = soln;
3943 503 : }
3944 :
3945 : void
3946 3541431 : NonlinearSystemBase::mortarConstraints(const Moose::ComputeType compute_type,
3947 : const std::set<TagID> & vector_tags,
3948 : const std::set<TagID> & matrix_tags)
3949 : {
3950 : parallel_object_only();
3951 :
3952 : try
3953 : {
3954 3549888 : for (auto & map_pr : _undisplaced_mortar_functors)
3955 8457 : map_pr.second(compute_type, vector_tags, matrix_tags);
3956 :
3957 3545696 : for (auto & map_pr : _displaced_mortar_functors)
3958 4265 : map_pr.second(compute_type, vector_tags, matrix_tags);
3959 : }
3960 0 : catch (MetaPhysicL::LogicError &)
3961 : {
3962 0 : mooseError(
3963 : "We caught a MetaPhysicL error in NonlinearSystemBase::mortarConstraints. This is very "
3964 : "likely due to AD not having a sufficiently large derivative container size. Please run "
3965 : "MOOSE configure with the '--with-derivative-size=<n>' option");
3966 0 : }
3967 3541431 : }
3968 :
3969 : void
3970 475 : NonlinearSystemBase::setupScalingData()
3971 : {
3972 475 : if (_auto_scaling_initd)
3973 0 : return;
3974 :
3975 : // Want the libMesh count of variables, not MOOSE, e.g. I don't care about array variable counts
3976 475 : const auto n_vars = system().n_vars();
3977 :
3978 475 : if (_scaling_group_variables.empty())
3979 : {
3980 464 : _var_to_group_var.reserve(n_vars);
3981 464 : _num_scaling_groups = n_vars;
3982 :
3983 1229 : for (const auto var_number : make_range(n_vars))
3984 765 : _var_to_group_var.emplace(var_number, var_number);
3985 : }
3986 : else
3987 : {
3988 11 : std::set<unsigned int> var_numbers, var_numbers_covered, var_numbers_not_covered;
3989 44 : for (const auto var_number : make_range(n_vars))
3990 33 : var_numbers.insert(var_number);
3991 :
3992 11 : _num_scaling_groups = _scaling_group_variables.size();
3993 :
3994 22 : for (const auto group_index : index_range(_scaling_group_variables))
3995 44 : for (const auto & var_name : _scaling_group_variables[group_index])
3996 : {
3997 33 : if (!hasVariable(var_name) && !hasScalarVariable(var_name))
3998 0 : mooseError("'",
3999 : var_name,
4000 : "', provided to the 'scaling_group_variables' parameter, does not exist in "
4001 : "the nonlinear system.");
4002 :
4003 : const MooseVariableBase & var =
4004 33 : hasVariable(var_name)
4005 22 : ? static_cast<MooseVariableBase &>(getVariable(0, var_name))
4006 55 : : static_cast<MooseVariableBase &>(getScalarVariable(0, var_name));
4007 33 : auto map_pair = _var_to_group_var.emplace(var.number(), group_index);
4008 33 : if (!map_pair.second)
4009 0 : mooseError("Variable ", var_name, " is contained in multiple scaling grouplings");
4010 33 : var_numbers_covered.insert(var.number());
4011 : }
4012 :
4013 11 : std::set_difference(var_numbers.begin(),
4014 : var_numbers.end(),
4015 : var_numbers_covered.begin(),
4016 : var_numbers_covered.end(),
4017 : std::inserter(var_numbers_not_covered, var_numbers_not_covered.begin()));
4018 :
4019 11 : _num_scaling_groups = _scaling_group_variables.size() + var_numbers_not_covered.size();
4020 :
4021 11 : auto index = static_cast<unsigned int>(_scaling_group_variables.size());
4022 11 : for (auto var_number : var_numbers_not_covered)
4023 0 : _var_to_group_var.emplace(var_number, index++);
4024 11 : }
4025 :
4026 475 : _variable_autoscaled.resize(n_vars, true);
4027 475 : const auto & number_to_var_map = _vars[0].numberToVariableMap();
4028 :
4029 475 : if (_ignore_variables_for_autoscaling.size())
4030 45 : for (const auto i : index_range(_variable_autoscaled))
4031 36 : if (std::find(_ignore_variables_for_autoscaling.begin(),
4032 : _ignore_variables_for_autoscaling.end(),
4033 72 : libmesh_map_find(number_to_var_map, i)->name()) !=
4034 72 : _ignore_variables_for_autoscaling.end())
4035 18 : _variable_autoscaled[i] = false;
4036 :
4037 475 : _auto_scaling_initd = true;
4038 : }
4039 :
4040 : bool
4041 1545 : NonlinearSystemBase::computeScaling()
4042 : {
4043 1545 : if (_compute_scaling_once && _computed_scaling)
4044 935 : return true;
4045 :
4046 610 : _console << "\nPerforming automatic scaling calculation\n" << std::endl;
4047 :
4048 3050 : TIME_SECTION("computeScaling", 3, "Computing Automatic Scaling");
4049 :
4050 : // It's funny but we need to assemble our vector of scaling factors here otherwise we will be
4051 : // applying scaling factors of 0 during Assembly of our scaling Jacobian
4052 610 : assembleScalingVector();
4053 :
4054 : // container for repeated access of element global dof indices
4055 610 : std::vector<dof_id_type> dof_indices;
4056 :
4057 610 : if (!_auto_scaling_initd)
4058 475 : setupScalingData();
4059 :
4060 610 : std::vector<Real> inverse_scaling_factors(_num_scaling_groups, 0);
4061 610 : std::vector<Real> resid_inverse_scaling_factors(_num_scaling_groups, 0);
4062 610 : std::vector<Real> jac_inverse_scaling_factors(_num_scaling_groups, 0);
4063 610 : auto & dof_map = dofMap();
4064 :
4065 : // what types of scaling do we want?
4066 610 : bool jac_scaling = _resid_vs_jac_scaling_param < 1. - TOLERANCE;
4067 610 : bool resid_scaling = _resid_vs_jac_scaling_param > TOLERANCE;
4068 :
4069 610 : const NumericVector<Number> & scaling_residual = RHS();
4070 :
4071 610 : if (jac_scaling)
4072 : {
4073 : // if (!_auto_scaling_initd)
4074 : // We need to reinit this when the number of dofs changes
4075 : // but there is no good way to track that
4076 : // In theory, it is the job of libmesh system to track this,
4077 : // but this special matrix is not owned by libMesh system
4078 : // Let us reinit eveytime since it is not expensive
4079 : {
4080 565 : auto init_vector = NumericVector<Number>::build(this->comm());
4081 565 : init_vector->init(system().n_dofs(), system().n_local_dofs(), /*fast=*/false, PARALLEL);
4082 :
4083 565 : _scaling_matrix->clear();
4084 565 : _scaling_matrix->init(*init_vector);
4085 565 : }
4086 :
4087 565 : _fe_problem.computingScalingJacobian(true);
4088 : // Dispatch to derived classes to ensure that we use the correct matrix tag
4089 565 : computeScalingJacobian();
4090 565 : _fe_problem.computingScalingJacobian(false);
4091 : }
4092 :
4093 610 : if (resid_scaling)
4094 : {
4095 45 : _fe_problem.computingScalingResidual(true);
4096 45 : _fe_problem.computingNonlinearResid(true);
4097 : // Dispatch to derived classes to ensure that we use the correct vector tag
4098 45 : computeScalingResidual();
4099 45 : _fe_problem.computingNonlinearResid(false);
4100 45 : _fe_problem.computingScalingResidual(false);
4101 : }
4102 :
4103 : // Did something bad happen during residual/Jacobian scaling computation?
4104 610 : if (_fe_problem.getFailNextNonlinearConvergenceCheck())
4105 0 : return false;
4106 :
4107 200736 : auto examine_dof_indices = [this,
4108 : jac_scaling,
4109 : resid_scaling,
4110 : &dof_map,
4111 : &jac_inverse_scaling_factors,
4112 : &resid_inverse_scaling_factors,
4113 : &scaling_residual](const auto & dof_indices, const auto var_number)
4114 : {
4115 1016826 : for (auto dof_index : dof_indices)
4116 816090 : if (dof_map.local_index(dof_index))
4117 : {
4118 803417 : if (jac_scaling)
4119 : {
4120 : // For now we will use the diagonal for determining scaling
4121 799904 : auto mat_value = (*_scaling_matrix)(dof_index, dof_index);
4122 799904 : auto & factor = jac_inverse_scaling_factors[_var_to_group_var[var_number]];
4123 799904 : factor = std::max(factor, std::abs(mat_value));
4124 : }
4125 803417 : if (resid_scaling)
4126 : {
4127 3513 : auto vec_value = scaling_residual(dof_index);
4128 3513 : auto & factor = resid_inverse_scaling_factors[_var_to_group_var[var_number]];
4129 3513 : factor = std::max(factor, std::abs(vec_value));
4130 : }
4131 : }
4132 200736 : };
4133 :
4134 : // Compute our scaling factors for the spatial field variables
4135 75644 : for (const auto & elem : _fe_problem.getCurrentAlgebraicElementRange())
4136 276564 : for (const auto i : make_range(system().n_vars()))
4137 201530 : if (_variable_autoscaled[i] && system().variable_type(i).family != SCALAR)
4138 : {
4139 200658 : dof_map.dof_indices(elem, dof_indices, i);
4140 200658 : examine_dof_indices(dof_indices, i);
4141 : }
4142 :
4143 1642 : for (const auto i : make_range(system().n_vars()))
4144 1032 : if (_variable_autoscaled[i] && system().variable_type(i).family == SCALAR)
4145 : {
4146 78 : dof_map.SCALAR_dof_indices(dof_indices, i);
4147 78 : examine_dof_indices(dof_indices, i);
4148 : }
4149 :
4150 610 : if (resid_scaling)
4151 45 : _communicator.max(resid_inverse_scaling_factors);
4152 610 : if (jac_scaling)
4153 565 : _communicator.max(jac_inverse_scaling_factors);
4154 :
4155 610 : if (jac_scaling && resid_scaling)
4156 0 : for (MooseIndex(inverse_scaling_factors) i = 0; i < inverse_scaling_factors.size(); ++i)
4157 : {
4158 : // Be careful not to take log(0)
4159 0 : if (!resid_inverse_scaling_factors[i])
4160 : {
4161 0 : if (!jac_inverse_scaling_factors[i])
4162 0 : inverse_scaling_factors[i] = 1;
4163 : else
4164 0 : inverse_scaling_factors[i] = jac_inverse_scaling_factors[i];
4165 : }
4166 0 : else if (!jac_inverse_scaling_factors[i])
4167 : // We know the resid is not zero
4168 0 : inverse_scaling_factors[i] = resid_inverse_scaling_factors[i];
4169 : else
4170 0 : inverse_scaling_factors[i] =
4171 0 : std::exp(_resid_vs_jac_scaling_param * std::log(resid_inverse_scaling_factors[i]) +
4172 0 : (1 - _resid_vs_jac_scaling_param) * std::log(jac_inverse_scaling_factors[i]));
4173 0 : }
4174 610 : else if (jac_scaling)
4175 565 : inverse_scaling_factors = jac_inverse_scaling_factors;
4176 45 : else if (resid_scaling)
4177 45 : inverse_scaling_factors = resid_inverse_scaling_factors;
4178 : else
4179 0 : mooseError("We shouldn't be calling this routine if we're not performing any scaling");
4180 :
4181 : // We have to make sure that our scaling values are not zero
4182 1620 : for (auto & scaling_factor : inverse_scaling_factors)
4183 1010 : if (scaling_factor == 0)
4184 36 : scaling_factor = 1;
4185 :
4186 : // Now flatten the group scaling factors to the individual variable scaling factors
4187 610 : std::vector<Real> flattened_inverse_scaling_factors(system().n_vars());
4188 1642 : for (const auto i : index_range(flattened_inverse_scaling_factors))
4189 1032 : flattened_inverse_scaling_factors[i] = inverse_scaling_factors[_var_to_group_var[i]];
4190 :
4191 : // Now set the scaling factors for the variables
4192 610 : applyScalingFactors(flattened_inverse_scaling_factors);
4193 610 : if (auto displaced_problem = _fe_problem.getDisplacedProblem().get())
4194 60 : displaced_problem->systemBaseNonlinear(number()).applyScalingFactors(
4195 : flattened_inverse_scaling_factors);
4196 :
4197 610 : _computed_scaling = true;
4198 610 : return true;
4199 610 : }
4200 :
4201 : void
4202 291417 : NonlinearSystemBase::assembleScalingVector()
4203 : {
4204 874251 : if (!hasVector("scaling_factors"))
4205 : // No variables have indicated they need scaling
4206 290873 : return;
4207 :
4208 1088 : auto & scaling_vector = getVector("scaling_factors");
4209 :
4210 544 : const auto & lm_mesh = _mesh.getMesh();
4211 544 : const auto & dof_map = dofMap();
4212 :
4213 544 : const auto & field_variables = _vars[0].fieldVariables();
4214 544 : const auto & scalar_variables = _vars[0].scalars();
4215 :
4216 544 : std::vector<dof_id_type> dof_indices;
4217 :
4218 544 : for (const Elem * const elem :
4219 12812 : as_range(lm_mesh.active_local_elements_begin(), lm_mesh.active_local_elements_end()))
4220 28956 : for (const auto * const field_var : field_variables)
4221 : {
4222 17232 : const auto & factors = field_var->arrayScalingFactor();
4223 34688 : for (const auto i : make_range(field_var->count()))
4224 : {
4225 17456 : dof_map.dof_indices(elem, dof_indices, field_var->number() + i);
4226 72152 : for (const auto dof : dof_indices)
4227 54696 : scaling_vector.set(dof, factors[i]);
4228 : }
4229 544 : }
4230 :
4231 646 : for (const auto * const scalar_var : scalar_variables)
4232 : {
4233 : mooseAssert(scalar_var->count() == 1,
4234 : "Scalar variables should always have only one component.");
4235 102 : dof_map.SCALAR_dof_indices(dof_indices, scalar_var->number());
4236 204 : for (const auto dof : dof_indices)
4237 102 : scaling_vector.set(dof, scalar_var->scalingFactor());
4238 : }
4239 :
4240 : // Parallel assemble
4241 544 : scaling_vector.close();
4242 :
4243 544 : if (auto * displaced_problem = _fe_problem.getDisplacedProblem().get())
4244 : // copy into the corresponding displaced system vector because they should be the exact same
4245 0 : displaced_problem->systemBaseNonlinear(number()).getVector("scaling_factors") = scaling_vector;
4246 544 : }
4247 :
4248 : bool
4249 290807 : NonlinearSystemBase::preSolve()
4250 : {
4251 : // Clear the iteration counters
4252 290807 : _current_l_its.clear();
4253 290807 : _current_nl_its = 0;
4254 :
4255 : // Initialize the solution vector using a predictor and known values from nodal bcs
4256 290807 : setInitialSolution();
4257 :
4258 : // Now that the initial solution has ben set, potentially perform a residual/Jacobian evaluation
4259 : // to determine variable scaling factors
4260 290807 : if (_automatic_scaling)
4261 : {
4262 1545 : const bool scaling_succeeded = computeScaling();
4263 1545 : if (!scaling_succeeded)
4264 0 : return false;
4265 : }
4266 :
4267 : // We do not know a priori what variable a global degree of freedom corresponds to, so we need a
4268 : // map from global dof to scaling factor. We just use a ghosted NumericVector for that mapping
4269 290807 : assembleScalingVector();
4270 :
4271 290807 : return true;
4272 : }
4273 :
4274 : void
4275 5058 : NonlinearSystemBase::destroyColoring()
4276 : {
4277 5058 : if (matrixFromColoring())
4278 108 : LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
4279 5058 : }
4280 :
4281 : FieldSplitPreconditionerBase &
4282 0 : NonlinearSystemBase::getFieldSplitPreconditioner()
4283 : {
4284 0 : if (!_fsp)
4285 0 : mooseError("No field split preconditioner is present for this system");
4286 :
4287 0 : return *_fsp;
4288 : }
|
