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