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