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 63676 : NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem,
115 : System & sys,
116 63676 : const std::string & name)
117 : : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
118 63676 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
119 63676 : _sys(sys),
120 63676 : _last_nl_rnorm(0.),
121 63676 : _current_nl_its(0),
122 63676 : _residual_ghosted(NULL),
123 63676 : _Re_time_tag(-1),
124 63676 : _Re_time(NULL),
125 63676 : _Re_non_time_tag(-1),
126 63676 : _Re_non_time(NULL),
127 63676 : _scalar_kernels(/*threaded=*/false),
128 63676 : _nodal_bcs(/*threaded=*/false),
129 63676 : _preset_nodal_bcs(/*threaded=*/false),
130 63676 : _ad_preset_nodal_bcs(/*threaded=*/false),
131 : #ifdef MOOSE_KOKKOS_ENABLED
132 42552 : _kokkos_kernels(/*threaded=*/false),
133 42552 : _kokkos_integrated_bcs(/*threaded=*/false),
134 42552 : _kokkos_nodal_bcs(/*threaded=*/false),
135 42552 : _kokkos_preset_nodal_bcs(/*threaded=*/false),
136 42552 : _kokkos_nodal_kernels(/*threaded=*/false),
137 : #endif
138 63676 : _general_dampers(/*threaded=*/false),
139 63676 : _splits(/*threaded=*/false),
140 63676 : _increment_vec(NULL),
141 63676 : _use_finite_differenced_preconditioner(false),
142 63676 : _fdcoloring(nullptr),
143 63676 : _fsp(nullptr),
144 63676 : _add_implicit_geometric_coupling_entries_to_jacobian(false),
145 63676 : _assemble_constraints_separately(false),
146 63676 : _need_residual_ghosted(false),
147 63676 : _debugging_residuals(false),
148 63676 : _doing_dg(false),
149 63676 : _n_iters(0),
150 63676 : _n_linear_iters(0),
151 63676 : _n_residual_evaluations(0),
152 63676 : _final_residual(0.),
153 63676 : _computing_pre_smo_residual(false),
154 63676 : _pre_smo_residual(0),
155 63676 : _initial_residual(0),
156 63676 : _use_pre_smo_residual(false),
157 63676 : _print_all_var_norms(false),
158 63676 : _has_save_in(false),
159 63676 : _has_diag_save_in(false),
160 63676 : _has_nodalbc_save_in(false),
161 63676 : _has_nodalbc_diag_save_in(false),
162 63676 : _computed_scaling(false),
163 63676 : _compute_scaling_once(true),
164 63676 : _resid_vs_jac_scaling_param(0),
165 63676 : _off_diagonals_in_auto_scaling(false),
166 509408 : _auto_scaling_initd(false)
167 : {
168 63676 : getResidualNonTimeVector();
169 : // Don't need to add the matrix - it already exists (for now)
170 63676 : _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 63676 : _fe_problem.addMatrixTag("TIME");
175 :
176 63676 : _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
177 :
178 63676 : _sys.identify_variable_groups(_fe_problem.identifyVariableGroupsInNL());
179 :
180 63676 : if (!_fe_problem.defaultGhosting())
181 : {
182 63594 : auto & dof_map = _sys.get_dof_map();
183 63594 : dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting());
184 63594 : dof_map.set_implicit_neighbor_dofs(false);
185 : }
186 63676 : }
187 :
188 59184 : NonlinearSystemBase::~NonlinearSystemBase() = default;
189 :
190 : void
191 61623 : NonlinearSystemBase::preInit()
192 : {
193 61623 : SolverSystem::preInit();
194 :
195 61623 : if (_fe_problem.hasDampers())
196 159 : setupDampers();
197 :
198 61623 : if (_residual_copy.get())
199 0 : _residual_copy->init(_sys.n_dofs(), false, SERIAL);
200 :
201 : #ifdef MOOSE_KOKKOS_ENABLED
202 41357 : if (_fe_problem.hasKokkosObjects())
203 741 : _sys.get_dof_map().full_sparsity_pattern_needed();
204 : #endif
205 61623 : }
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 60693 : NonlinearSystemBase::initialSetup()
227 : {
228 303465 : TIME_SECTION("nlInitialSetup", 2, "Setting Up Nonlinear System");
229 :
230 60693 : SolverSystem::initialSetup();
231 :
232 : {
233 303465 : TIME_SECTION("kernelsInitialSetup", 2, "Setting Up Kernels/BCs/Constraints");
234 :
235 127047 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
236 : {
237 66354 : _kernels.initialSetup(tid);
238 66354 : _nodal_kernels.initialSetup(tid);
239 66354 : _dirac_kernels.initialSetup(tid);
240 66354 : if (_doing_dg)
241 1225 : _dg_kernels.initialSetup(tid);
242 66354 : _interface_kernels.initialSetup(tid);
243 :
244 66354 : _element_dampers.initialSetup(tid);
245 66354 : _nodal_dampers.initialSetup(tid);
246 66354 : _integrated_bcs.initialSetup(tid);
247 :
248 66354 : if (_fe_problem.haveFV())
249 : {
250 4338 : std::vector<FVElementalKernel *> fv_elemental_kernels;
251 4338 : _fe_problem.theWarehouse()
252 8676 : .query()
253 4338 : .template condition<AttribSystem>("FVElementalKernel")
254 4338 : .template condition<AttribThread>(tid)
255 4338 : .queryInto(fv_elemental_kernels);
256 :
257 8340 : for (auto * fv_kernel : fv_elemental_kernels)
258 4002 : fv_kernel->initialSetup();
259 :
260 4338 : std::vector<FVFluxKernel *> fv_flux_kernels;
261 4338 : _fe_problem.theWarehouse()
262 8676 : .query()
263 4338 : .template condition<AttribSystem>("FVFluxKernel")
264 4338 : .template condition<AttribThread>(tid)
265 4338 : .queryInto(fv_flux_kernels);
266 :
267 9795 : for (auto * fv_kernel : fv_flux_kernels)
268 5457 : fv_kernel->initialSetup();
269 4338 : }
270 : }
271 :
272 60693 : _scalar_kernels.initialSetup();
273 60693 : _constraints.initialSetup();
274 60693 : _general_dampers.initialSetup();
275 60693 : _nodal_bcs.initialSetup();
276 60693 : _preset_nodal_bcs.residualSetup();
277 60693 : _ad_preset_nodal_bcs.residualSetup();
278 :
279 : #ifdef MOOSE_KOKKOS_ENABLED
280 40872 : _kokkos_kernels.initialSetup();
281 40872 : _kokkos_nodal_kernels.initialSetup();
282 40872 : _kokkos_integrated_bcs.initialSetup();
283 40872 : _kokkos_nodal_bcs.initialSetup();
284 : #endif
285 60693 : }
286 :
287 : {
288 303465 : TIME_SECTION("mortarSetup", 2, "Initializing Mortar Interfaces");
289 :
290 121386 : auto create_mortar_functors = [this](const bool displaced)
291 : {
292 : // go over mortar interfaces and construct functors
293 121386 : const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
294 122406 : for (const auto & [primary_secondary_boundary_pair, mortar_generation_ptr] :
295 243792 : 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 121386 : };
319 :
320 60693 : create_mortar_functors(false);
321 60693 : create_mortar_functors(true);
322 60693 : }
323 :
324 60693 : 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 60693 : if (_preconditioner)
333 14709 : _preconditioner->initialSetup();
334 60693 : }
335 :
336 : void
337 291885 : NonlinearSystemBase::timestepSetup()
338 : {
339 291885 : SolverSystem::timestepSetup();
340 :
341 609724 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
342 : {
343 317839 : _kernels.timestepSetup(tid);
344 317839 : _nodal_kernels.timestepSetup(tid);
345 317839 : _dirac_kernels.timestepSetup(tid);
346 317839 : if (_doing_dg)
347 1799 : _dg_kernels.timestepSetup(tid);
348 317839 : _interface_kernels.timestepSetup(tid);
349 317839 : _element_dampers.timestepSetup(tid);
350 317839 : _nodal_dampers.timestepSetup(tid);
351 317839 : _integrated_bcs.timestepSetup(tid);
352 :
353 317839 : if (_fe_problem.haveFV())
354 : {
355 17448 : std::vector<FVFluxBC *> bcs;
356 17448 : _fe_problem.theWarehouse()
357 34896 : .query()
358 17448 : .template condition<AttribSystem>("FVFluxBC")
359 17448 : .template condition<AttribThread>(tid)
360 17448 : .queryInto(bcs);
361 :
362 17448 : std::vector<FVInterfaceKernel *> iks;
363 17448 : _fe_problem.theWarehouse()
364 34896 : .query()
365 17448 : .template condition<AttribSystem>("FVInterfaceKernel")
366 17448 : .template condition<AttribThread>(tid)
367 17448 : .queryInto(iks);
368 :
369 17448 : std::vector<FVFluxKernel *> kernels;
370 17448 : _fe_problem.theWarehouse()
371 34896 : .query()
372 17448 : .template condition<AttribSystem>("FVFluxKernel")
373 17448 : .template condition<AttribThread>(tid)
374 17448 : .queryInto(kernels);
375 :
376 34465 : for (auto * bc : bcs)
377 17017 : bc->timestepSetup();
378 17846 : for (auto * ik : iks)
379 398 : ik->timestepSetup();
380 38624 : for (auto * kernel : kernels)
381 21176 : kernel->timestepSetup();
382 17448 : }
383 : }
384 291885 : _scalar_kernels.timestepSetup();
385 291885 : _constraints.timestepSetup();
386 291885 : _general_dampers.timestepSetup();
387 291885 : _nodal_bcs.timestepSetup();
388 291885 : _preset_nodal_bcs.timestepSetup();
389 291885 : _ad_preset_nodal_bcs.timestepSetup();
390 :
391 : #ifdef MOOSE_KOKKOS_ENABLED
392 192369 : _kokkos_kernels.timestepSetup();
393 192369 : _kokkos_nodal_kernels.timestepSetup();
394 192369 : _kokkos_integrated_bcs.timestepSetup();
395 192369 : _kokkos_nodal_bcs.timestepSetup();
396 : #endif
397 291885 : }
398 :
399 : void
400 1951317 : NonlinearSystemBase::customSetup(const ExecFlagType & exec_type)
401 : {
402 1951317 : SolverSystem::customSetup(exec_type);
403 :
404 4077703 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
405 : {
406 2126386 : _kernels.customSetup(exec_type, tid);
407 2126386 : _nodal_kernels.customSetup(exec_type, tid);
408 2126386 : _dirac_kernels.customSetup(exec_type, tid);
409 2126386 : if (_doing_dg)
410 12322 : _dg_kernels.customSetup(exec_type, tid);
411 2126386 : _interface_kernels.customSetup(exec_type, tid);
412 2126386 : _element_dampers.customSetup(exec_type, tid);
413 2126386 : _nodal_dampers.customSetup(exec_type, tid);
414 2126386 : _integrated_bcs.customSetup(exec_type, tid);
415 :
416 2126386 : if (_fe_problem.haveFV())
417 : {
418 128935 : std::vector<FVFluxBC *> bcs;
419 128935 : _fe_problem.theWarehouse()
420 257870 : .query()
421 128935 : .template condition<AttribSystem>("FVFluxBC")
422 128935 : .template condition<AttribThread>(tid)
423 128935 : .queryInto(bcs);
424 :
425 128935 : std::vector<FVInterfaceKernel *> iks;
426 128935 : _fe_problem.theWarehouse()
427 257870 : .query()
428 128935 : .template condition<AttribSystem>("FVInterfaceKernel")
429 128935 : .template condition<AttribThread>(tid)
430 128935 : .queryInto(iks);
431 :
432 128935 : std::vector<FVFluxKernel *> kernels;
433 128935 : _fe_problem.theWarehouse()
434 257870 : .query()
435 128935 : .template condition<AttribSystem>("FVFluxKernel")
436 128935 : .template condition<AttribThread>(tid)
437 128935 : .queryInto(kernels);
438 :
439 232415 : for (auto * bc : bcs)
440 103480 : bc->customSetup(exec_type);
441 173154 : for (auto * ik : iks)
442 44219 : ik->customSetup(exec_type);
443 306653 : for (auto * kernel : kernels)
444 177718 : kernel->customSetup(exec_type);
445 128935 : }
446 : }
447 1951317 : _scalar_kernels.customSetup(exec_type);
448 1951317 : _constraints.customSetup(exec_type);
449 1951317 : _general_dampers.customSetup(exec_type);
450 1951317 : _nodal_bcs.customSetup(exec_type);
451 1951317 : _preset_nodal_bcs.customSetup(exec_type);
452 1951317 : _ad_preset_nodal_bcs.customSetup(exec_type);
453 :
454 : #ifdef MOOSE_KOKKOS_ENABLED
455 1285328 : _kokkos_kernels.customSetup(exec_type);
456 1285328 : _kokkos_nodal_kernels.customSetup(exec_type);
457 1285328 : _kokkos_integrated_bcs.customSetup(exec_type);
458 1285328 : _kokkos_nodal_bcs.customSetup(exec_type);
459 : #endif
460 1951317 : }
461 :
462 : void
463 348072 : NonlinearSystemBase::setupDM()
464 : {
465 348072 : if (_fsp)
466 129 : _fsp->setupDM();
467 348072 : }
468 :
469 : void
470 83363 : NonlinearSystemBase::addKernel(const std::string & kernel_name,
471 : const std::string & name,
472 : InputParameters & parameters)
473 : {
474 174277 : 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 91118 : _factory.create<KernelBase>(kernel_name, name, parameters, tid);
479 90938 : _kernels.addObject(kernel, tid);
480 90918 : postAddResidualObject(*kernel);
481 : // Add to theWarehouse, a centralized storage for all moose objects
482 90914 : _fe_problem.theWarehouse().add(kernel);
483 90914 : }
484 :
485 83159 : if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
486 130 : _has_save_in = true;
487 83159 : if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
488 104 : _has_diag_save_in = true;
489 83159 : }
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 915 : NonlinearSystemBase::addDiracKernel(const std::string & kernel_name,
648 : const std::string & name,
649 : InputParameters & parameters)
650 : {
651 1906 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
652 : {
653 : std::shared_ptr<DiracKernelBase> kernel =
654 999 : _factory.create<DiracKernelBase>(kernel_name, name, parameters, tid);
655 991 : postAddResidualObject(*kernel);
656 991 : _dirac_kernels.addObject(kernel, tid);
657 : // Add to theWarehouse, a centralized storage for all moose objects
658 991 : _fe_problem.theWarehouse().add(kernel);
659 991 : }
660 907 : }
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 316761 : NonlinearSystemBase::shouldEvaluatePreSMOResidual() const
753 : {
754 316761 : 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 305954 : if (_app.parameters().get<bool>("use_legacy_initial_residual_evaluation_behavior"))
764 0 : return true;
765 :
766 305954 : return _use_pre_smo_residual;
767 : }
768 :
769 : Real
770 394915 : NonlinearSystemBase::referenceResidual() const
771 : {
772 394915 : 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 395027 : NonlinearSystemBase::initialResidual() const
786 : {
787 395027 : return _initial_residual;
788 : }
789 :
790 : void
791 313938 : NonlinearSystemBase::setInitialResidual(Real r)
792 : {
793 313938 : _initial_residual = r;
794 313938 : }
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 3339314 : NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
830 : {
831 : parallel_object_only();
832 :
833 10017942 : TIME_SECTION("nl::computeResidualTags", 5);
834 :
835 3339314 : _fe_problem.setCurrentNonlinearSystem(number());
836 3339314 : _fe_problem.setCurrentlyComputingResidual(true);
837 :
838 3339314 : bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
839 :
840 3339314 : _n_residual_evaluations++;
841 :
842 : // not suppose to do anythin on matrix
843 3339314 : deactivateAllMatrixTags();
844 :
845 3339314 : FloatingPointExceptionGuard fpe_guard(_app);
846 :
847 3339314 : 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 3339314 : zeroTaggedVectors(tags);
858 3339314 : computeResidualInternal(tags);
859 3339149 : closeTaggedVectors(tags);
860 :
861 3339149 : if (required_residual)
862 : {
863 3307237 : auto & residual = getVector(residualVectorTag());
864 3307237 : if (!_time_integrators.empty())
865 : {
866 5817365 : for (auto & ti : _time_integrators)
867 2917108 : ti->postResidual(residual);
868 : }
869 : else
870 406980 : residual += *_Re_non_time;
871 3307237 : residual.close();
872 : }
873 3339149 : if (_fe_problem.computingScalingResidual())
874 : // We don't want to do nodal bcs or anything else
875 50 : return;
876 :
877 3339099 : computeNodalBCs(tags);
878 3339099 : closeTaggedVectors(tags);
879 :
880 : // If we are debugging residuals we need one more assignment to have the ghosted copy up to
881 : // date
882 3339099 : 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 3339099 : if (_has_nodalbc_save_in)
891 171 : _fe_problem.getAuxiliarySystem().solution().close();
892 3339099 : 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 3339220 : activateAllMatrixTags();
904 :
905 3339220 : _fe_problem.setCurrentlyComputingResidual(false);
906 3339320 : }
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 262599 : NonlinearSystemBase::onTimestepBegin()
949 : {
950 525916 : for (auto & ti : _time_integrators)
951 263317 : ti->preSolve();
952 262599 : if (_predictor.get())
953 387 : _predictor->timestepSetup();
954 262599 : }
955 :
956 : void
957 317896 : NonlinearSystemBase::setInitialSolution()
958 : {
959 317896 : deactivateAllMatrixTags();
960 :
961 317896 : NumericVector<Number> & initial_solution(solution());
962 317896 : 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 1589480 : TIME_SECTION("initialBCs", 2, "Applying BCs To Initial Condition");
978 :
979 317896 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
980 15570900 : for (const auto & bnode : bnd_nodes)
981 : {
982 15253004 : BoundaryID boundary_id = bnode->_bnd_id;
983 15253004 : Node * node = bnode->_node;
984 :
985 15253004 : if (node->processor_id() == processor_id())
986 : {
987 11703838 : bool has_preset_nodal_bcs = _preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
988 11703838 : bool has_ad_preset_nodal_bcs = _ad_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
989 :
990 : // reinit variables in nodes
991 11703838 : if (has_preset_nodal_bcs || has_ad_preset_nodal_bcs)
992 3685445 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
993 :
994 11703838 : if (has_preset_nodal_bcs)
995 : {
996 3648404 : const auto & preset_bcs = _preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
997 7955880 : for (const auto & preset_bc : preset_bcs)
998 4307476 : preset_bc->computeValue(initial_solution);
999 : }
1000 11703838 : 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 317896 : }
1009 :
1010 : #ifdef MOOSE_KOKKOS_ENABLED
1011 209372 : if (_kokkos_preset_nodal_bcs.hasObjects())
1012 1889 : setKokkosInitialSolution();
1013 : #endif
1014 :
1015 317896 : _sys.solution->close();
1016 317896 : update();
1017 :
1018 : // Set constraint secondary values
1019 317896 : setConstraintSecondaryValues(initial_solution, false);
1020 :
1021 317896 : if (_fe_problem.getDisplacedProblem())
1022 35761 : setConstraintSecondaryValues(initial_solution, true);
1023 317896 : }
1024 :
1025 : void
1026 121 : NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
1027 : {
1028 121 : _predictor = predictor;
1029 121 : }
1030 :
1031 : void
1032 6216003 : NonlinearSystemBase::subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
1033 : {
1034 6216003 : SolverSystem::subdomainSetup();
1035 :
1036 6216003 : _kernels.subdomainSetup(subdomain, tid);
1037 6216003 : _nodal_kernels.subdomainSetup(subdomain, tid);
1038 6216003 : _element_dampers.subdomainSetup(subdomain, tid);
1039 6216003 : _nodal_dampers.subdomainSetup(subdomain, tid);
1040 6216003 : }
1041 :
1042 : NumericVector<Number> &
1043 32071 : NonlinearSystemBase::getResidualTimeVector()
1044 : {
1045 32071 : if (!_Re_time)
1046 : {
1047 31965 : _Re_time_tag = _fe_problem.addVectorTag("TIME");
1048 :
1049 : // Most applications don't need the expense of ghosting
1050 31965 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1051 31965 : _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 32071 : return *_Re_time;
1063 : }
1064 :
1065 : NumericVector<Number> &
1066 95954 : NonlinearSystemBase::getResidualNonTimeVector()
1067 : {
1068 95954 : if (!_Re_non_time)
1069 : {
1070 63676 : _Re_non_time_tag = _fe_problem.addVectorTag("NONTIME");
1071 :
1072 : // Most applications don't need the expense of ghosting
1073 63676 : ParallelType ptype = _need_residual_ghosted ? GHOSTED : PARALLEL;
1074 63676 : _Re_non_time = &addVector(_Re_non_time_tag, false, ptype);
1075 : }
1076 32278 : 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 95954 : 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 23409 : NonlinearSystemBase::enforceNodalConstraintsResidual(NumericVector<Number> & residual)
1107 : {
1108 23409 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1109 23409 : residual.close();
1110 23409 : 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 23409 : }
1129 :
1130 : void
1131 3667 : NonlinearSystemBase::enforceNodalConstraintsJacobian()
1132 : {
1133 3667 : if (!hasMatrix(systemMatrixTag()))
1134 0 : mooseError(" A system matrix is required");
1135 :
1136 3667 : auto & jacobian = getMatrix(systemMatrixTag());
1137 3667 : THREAD_ID tid = 0; // constraints are going to be done single-threaded
1138 :
1139 3667 : 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 3667 : }
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 353657 : 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 389418 : : static_cast<SubProblem &>(_fe_problem);
1230 353657 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1231 :
1232 353657 : bool constraints_applied = false;
1233 :
1234 413926 : 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 353657 : std::set<dof_id_type> unique_secondary_node_ids;
1304 :
1305 793733 : for (const auto & secondary_id : _mesh.meshSubdomains())
1306 : {
1307 1174610 : for (const auto & primary_id : _mesh.meshSubdomains())
1308 : {
1309 734534 : 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 353657 : _communicator.max(constraints_applied);
1354 :
1355 353657 : if (constraints_applied)
1356 : {
1357 719 : solution.close();
1358 719 : update();
1359 : }
1360 353657 : }
1361 :
1362 : void
1363 28287 : NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool displaced)
1364 : {
1365 : // Make sure the residual is in a good state
1366 28287 : 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 33165 : : static_cast<SubProblem &>(_fe_problem);
1374 28287 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1375 :
1376 : bool constraints_applied;
1377 28287 : bool residual_has_inserted_values = false;
1378 28287 : if (!_assemble_constraints_separately)
1379 28287 : constraints_applied = false;
1380 52353 : 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 28287 : if (!_assemble_constraints_separately)
1509 : {
1510 28287 : _communicator.max(constraints_applied);
1511 :
1512 28287 : 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 28287 : THREAD_ID tid = 0;
1530 28287 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
1531 28287 : 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 28287 : std::set<dof_id_type> unique_secondary_node_ids;
1577 :
1578 28287 : constraints_applied = false;
1579 28287 : residual_has_inserted_values = false;
1580 28287 : bool has_writable_variables = false;
1581 97398 : for (const auto & secondary_id : _mesh.meshSubdomains())
1582 : {
1583 296138 : for (const auto & primary_id : _mesh.meshSubdomains())
1584 : {
1585 227027 : 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 28287 : _communicator.max(constraints_applied);
1647 :
1648 28287 : 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 28287 : _communicator.max(has_writable_variables);
1663 :
1664 28287 : 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 28287 : _fe_problem.addCachedResidual(0);
1676 28287 : }
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 3435321 : NonlinearSystemBase::residualSetup()
1729 : {
1730 10305963 : TIME_SECTION("residualSetup", 3);
1731 :
1732 3435321 : SolverSystem::residualSetup();
1733 :
1734 7178550 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1735 : {
1736 3743229 : _kernels.residualSetup(tid);
1737 3743229 : _nodal_kernels.residualSetup(tid);
1738 3743229 : _dirac_kernels.residualSetup(tid);
1739 3743229 : if (_doing_dg)
1740 38033 : _dg_kernels.residualSetup(tid);
1741 3743229 : _interface_kernels.residualSetup(tid);
1742 3743229 : _element_dampers.residualSetup(tid);
1743 3743229 : _nodal_dampers.residualSetup(tid);
1744 3743229 : _integrated_bcs.residualSetup(tid);
1745 : }
1746 3435321 : _scalar_kernels.residualSetup();
1747 3435321 : _constraints.residualSetup();
1748 3435321 : _general_dampers.residualSetup();
1749 3435321 : _nodal_bcs.residualSetup();
1750 3435321 : _preset_nodal_bcs.residualSetup();
1751 3435321 : _ad_preset_nodal_bcs.residualSetup();
1752 :
1753 : #ifdef MOOSE_KOKKOS_ENABLED
1754 2255785 : _kokkos_kernels.residualSetup();
1755 2255785 : _kokkos_nodal_kernels.residualSetup();
1756 2255785 : _kokkos_integrated_bcs.residualSetup();
1757 2255785 : _kokkos_nodal_bcs.residualSetup();
1758 : #endif
1759 :
1760 : // Avoid recursion
1761 3435321 : if (this == &_fe_problem.currentNonlinearSystem())
1762 3342781 : _fe_problem.residualSetup();
1763 3435321 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
1764 3435321 : }
1765 :
1766 : void
1767 3339314 : NonlinearSystemBase::computeResidualInternal(const std::set<TagID> & tags)
1768 : {
1769 : parallel_object_only();
1770 :
1771 10017942 : TIME_SECTION("computeResidualInternal", 3);
1772 :
1773 3339314 : residualSetup();
1774 :
1775 : #ifdef MOOSE_KOKKOS_ENABLED
1776 2191652 : if (_fe_problem.hasKokkosResidualObjects())
1777 26756 : computeKokkosResidual(tags);
1778 : #endif
1779 :
1780 3339314 : 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 3339314 : std::vector<UserObject *> uos;
1785 3339314 : _fe_problem.theWarehouse()
1786 6678628 : .query()
1787 3339314 : .condition<AttribSystem>("UserObject")
1788 3339314 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1789 3339314 : .queryInto(uos);
1790 3339314 : for (auto & uo : uos)
1791 0 : uo->residualSetup();
1792 3339314 : for (auto & uo : uos)
1793 : {
1794 0 : uo->initialize();
1795 0 : uo->execute();
1796 0 : uo->finalize();
1797 : }
1798 :
1799 : // reinit scalar variables
1800 6978317 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1801 3639003 : _fe_problem.reinitScalars(tid);
1802 :
1803 : // residual contributions from the domain
1804 : PARALLEL_TRY
1805 : {
1806 10017942 : TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1807 :
1808 3339314 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
1809 :
1810 3339314 : ComputeResidualThread cr(_fe_problem, tags);
1811 3339314 : 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 3339290 : if (_fe_problem.haveFV())
1819 : {
1820 : ComputeFVFluxResidualThread<FVRange> fvr(
1821 150860 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
1822 150860 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
1823 150860 : Threads::parallel_reduce(faces, fvr);
1824 150860 : }
1825 3339282 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
1826 3339282 : 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 3339282 : }
1834 :
1835 3339282 : unsigned int n_threads = libMesh::n_threads();
1836 6978245 : for (unsigned int i = 0; i < n_threads;
1837 : i++) // Add any cached residuals that might be hanging around
1838 3638963 : _fe_problem.addCachedResidual(i);
1839 3339282 : }
1840 3339282 : 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 3339161 : 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 3339161 : PARALLEL_CATCH;
1887 :
1888 : // residual contributions from Block NodalKernels
1889 : PARALLEL_TRY
1890 : {
1891 3339161 : 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 3339161 : PARALLEL_CATCH;
1913 :
1914 3339161 : 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 3339111 : 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 3339111 : PARALLEL_CATCH;
1940 :
1941 3339111 : mortarConstraints(Moose::ComputeType::Residual, tags, {});
1942 :
1943 3339111 : if (_residual_copy.get())
1944 : {
1945 0 : _Re_non_time->close();
1946 0 : _Re_non_time->localize(*_residual_copy);
1947 : }
1948 :
1949 3339111 : 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 3339111 : PARALLEL_TRY { computeDiracContributions(tags, false); }
1957 3339099 : PARALLEL_CATCH;
1958 :
1959 3339099 : if (_fe_problem._has_constraints)
1960 : {
1961 23409 : PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1962 23409 : PARALLEL_CATCH;
1963 23409 : _Re_non_time->close();
1964 : }
1965 :
1966 : // Add in Residual contributions from other Constraints
1967 3339099 : if (_fe_problem._has_constraints)
1968 : {
1969 : PARALLEL_TRY
1970 : {
1971 : // Undisplaced Constraints
1972 23409 : constraintResiduals(*_Re_non_time, false);
1973 :
1974 : // Displaced Constraints
1975 23409 : if (_fe_problem.getDisplacedProblem())
1976 4878 : constraintResiduals(*_Re_non_time, true);
1977 :
1978 23409 : if (_fe_problem.computingNonlinearResid())
1979 10705 : _constraints.residualEnd();
1980 : }
1981 23409 : PARALLEL_CATCH;
1982 23409 : _Re_non_time->close();
1983 : }
1984 :
1985 : // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
1986 : // counters
1987 3339099 : _app.solutionInvalidity().syncIteration();
1988 3339099 : _app.solutionInvalidity().solutionInvalidAccumulation();
1989 3339612 : }
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 3339099 : 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 3339099 : if (_has_save_in)
2116 309 : _fe_problem.getAuxiliarySystem().solution().close();
2117 :
2118 : // Select nodal kernels
2119 : MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2120 :
2121 3339099 : if (tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL) || !tags.size())
2122 3300218 : 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 3339099 : if (!nbc_warehouse->size())
2130 450654 : return;
2131 :
2132 : PARALLEL_TRY
2133 : {
2134 2888445 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
2135 :
2136 2888445 : if (!bnd_nodes.empty())
2137 : {
2138 8664402 : TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2139 :
2140 151131426 : for (const auto & bnode : bnd_nodes)
2141 : {
2142 148243292 : BoundaryID boundary_id = bnode->_bnd_id;
2143 148243292 : Node * node = bnode->_node;
2144 :
2145 262271804 : if (node->processor_id() == processor_id() &&
2146 114028512 : nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
2147 : {
2148 : // reinit variables in nodes
2149 58384845 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
2150 :
2151 58384845 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2152 122399114 : for (const auto & nbc : bcs)
2153 64014269 : if (nbc->shouldApply())
2154 64012865 : nbc->computeResidual();
2155 : }
2156 : }
2157 2888134 : }
2158 : }
2159 2888445 : PARALLEL_CATCH;
2160 :
2161 2888445 : if (_Re_time)
2162 2565727 : _Re_time->close();
2163 2888445 : _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 4763 : NonlinearSystemBase::constraintJacobians(const SparseMatrix<Number> & jacobian_to_view,
2351 : bool displaced)
2352 : {
2353 4763 : if (!hasMatrix(systemMatrixTag()))
2354 0 : mooseError("A system matrix is required");
2355 :
2356 4763 : auto & jacobian = getMatrix(systemMatrixTag());
2357 :
2358 4763 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2359 0 : LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2360 : MAT_NEW_NONZERO_ALLOCATION_ERR,
2361 : PETSC_FALSE));
2362 4763 : if (_fe_problem.ignoreZerosInJacobian())
2363 0 : LibmeshPetscCall(MatSetOption(
2364 : static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2365 :
2366 4763 : 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 5859 : : static_cast<SubProblem &>(_fe_problem);
2374 4763 : const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
2375 :
2376 : bool constraints_applied;
2377 4763 : if (!_assemble_constraints_separately)
2378 4763 : constraints_applied = false;
2379 6829 : 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 4763 : if (!_assemble_constraints_separately)
2558 : {
2559 : // See if constraints were applied anywhere
2560 4763 : _communicator.max(constraints_applied);
2561 :
2562 4763 : 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 4763 : THREAD_ID tid = 0;
2577 : // go over element-element constraint interface
2578 4763 : const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
2579 4763 : 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 4763 : std::set<dof_id_type> unique_secondary_node_ids;
2625 4763 : constraints_applied = false;
2626 20461 : for (const auto & secondary_id : _mesh.meshSubdomains())
2627 : {
2628 77282 : for (const auto & primary_id : _mesh.meshSubdomains())
2629 : {
2630 61584 : 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 4763 : _communicator.max(constraints_applied);
2745 :
2746 4763 : 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 4763 : }
2759 :
2760 : void
2761 527890 : NonlinearSystemBase::computeScalarKernelsJacobians(const std::set<TagID> & tags)
2762 : {
2763 : MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
2764 :
2765 527890 : if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2766 518066 : 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 527890 : 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 527890 : }
2808 :
2809 : void
2810 543950 : NonlinearSystemBase::jacobianSetup()
2811 : {
2812 543950 : SolverSystem::jacobianSetup();
2813 :
2814 1138082 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2815 : {
2816 594132 : _kernels.jacobianSetup(tid);
2817 594132 : _nodal_kernels.jacobianSetup(tid);
2818 594132 : _dirac_kernels.jacobianSetup(tid);
2819 594132 : if (_doing_dg)
2820 2998 : _dg_kernels.jacobianSetup(tid);
2821 594132 : _interface_kernels.jacobianSetup(tid);
2822 594132 : _element_dampers.jacobianSetup(tid);
2823 594132 : _nodal_dampers.jacobianSetup(tid);
2824 594132 : _integrated_bcs.jacobianSetup(tid);
2825 : }
2826 543950 : _scalar_kernels.jacobianSetup();
2827 543950 : _constraints.jacobianSetup();
2828 543950 : _general_dampers.jacobianSetup();
2829 543950 : _nodal_bcs.jacobianSetup();
2830 543950 : _preset_nodal_bcs.jacobianSetup();
2831 543950 : _ad_preset_nodal_bcs.jacobianSetup();
2832 :
2833 : #ifdef MOOSE_KOKKOS_ENABLED
2834 355891 : _kokkos_kernels.jacobianSetup();
2835 355891 : _kokkos_nodal_kernels.jacobianSetup();
2836 355891 : _kokkos_integrated_bcs.jacobianSetup();
2837 355891 : _kokkos_nodal_bcs.jacobianSetup();
2838 : #endif
2839 :
2840 : // Avoid recursion
2841 543950 : if (this == &_fe_problem.currentNonlinearSystem())
2842 527890 : _fe_problem.jacobianSetup();
2843 543950 : _app.solutionInvalidity().resetSolutionInvalidCurrentIteration();
2844 543950 : }
2845 :
2846 : void
2847 527890 : NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
2848 : {
2849 1583670 : TIME_SECTION("computeJacobianInternal", 3);
2850 :
2851 527890 : _fe_problem.setCurrentNonlinearSystem(number());
2852 :
2853 : // Make matrix ready to use
2854 527890 : activateAllMatrixTags();
2855 :
2856 1575832 : for (auto tag : tags)
2857 : {
2858 1047942 : if (!hasMatrix(tag))
2859 518066 : continue;
2860 :
2861 529876 : auto & jacobian = getMatrix(tag);
2862 : // Necessary for speed
2863 529876 : if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2864 : {
2865 528590 : LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2866 : MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2867 : PETSC_TRUE));
2868 528590 : if (!_fe_problem.errorOnJacobianNonzeroReallocation())
2869 13371 : LibmeshPetscCall(
2870 : MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2871 528590 : 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 527890 : jacobianSetup();
2879 :
2880 : #ifdef MOOSE_KOKKOS_ENABLED
2881 345197 : 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 527890 : std::vector<UserObject *> uos;
2888 527890 : _fe_problem.theWarehouse()
2889 1055780 : .query()
2890 527890 : .condition<AttribSystem>("UserObject")
2891 527890 : .condition<AttribExecOns>(EXEC_PRE_KERNELS)
2892 527890 : .queryInto(uos);
2893 527890 : for (auto & uo : uos)
2894 0 : uo->jacobianSetup();
2895 527890 : for (auto & uo : uos)
2896 : {
2897 0 : uo->initialize();
2898 0 : uo->execute();
2899 0 : uo->finalize();
2900 : }
2901 :
2902 : // reinit scalar variables
2903 1104589 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2904 576699 : _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 527890 : computeScalarKernelsJacobians(tags);
2912 :
2913 : // Block restricted Nodal Kernels
2914 527890 : 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 527890 : 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 37610 : _fe_problem, this->number(), tags, /*on_displaced=*/false);
2933 37610 : FVRange faces(_fe_problem.mesh().ownedFaceInfoBegin(), _fe_problem.mesh().ownedFaceInfoEnd());
2934 37610 : Threads::parallel_reduce(faces, fvj);
2935 37610 : }
2936 527890 : if (auto displaced_problem = _fe_problem.getDisplacedProblem();
2937 527890 : 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 527890 : }
2945 :
2946 527890 : mortarConstraints(Moose::ComputeType::Jacobian, {}, tags);
2947 :
2948 : // Get our element range for looping over
2949 527890 : const ConstElemRange & elem_range = _fe_problem.getCurrentAlgebraicElementRange();
2950 :
2951 527890 : 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 527288 : 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 93867 : default:
3004 : case Moose::COUPLING_CUSTOM:
3005 : {
3006 93867 : ComputeFullJacobianThread cj(_fe_problem, tags);
3007 93867 : Threads::parallel_reduce(elem_range, cj);
3008 93867 : unsigned int n_threads = libMesh::n_threads();
3009 :
3010 194455 : for (unsigned int i = 0; i < n_threads; i++)
3011 100592 : _fe_problem.addCachedJacobian(i);
3012 :
3013 : // Boundary restricted Nodal Kernels
3014 93863 : 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 93867 : }
3026 93863 : break;
3027 : }
3028 :
3029 527276 : 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 1053928 : if ((_fe_problem.restoreOriginalNonzeroPattern() || first) &&
3035 526656 : _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 527272 : 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 527170 : 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 3667 : auto & system_matrix = getMatrix(systemMatrixTag());
3056 : #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
3057 : SparseMatrix<Number> * view_jac_ptr;
3058 3667 : std::unique_ptr<SparseMatrix<Number>> hash_copy;
3059 3667 : 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 1358 : view_jac_ptr = &system_matrix;
3066 3667 : auto & jacobian_to_view = *view_jac_ptr;
3067 : #else
3068 : auto & jacobian_to_view = system_matrix;
3069 : #endif
3070 3667 : if (&jacobian_to_view == &system_matrix)
3071 1358 : system_matrix.close();
3072 :
3073 : // Nodal Constraints
3074 3667 : enforceNodalConstraintsJacobian();
3075 :
3076 : // Undisplaced Constraints
3077 3667 : constraintJacobians(jacobian_to_view, false);
3078 :
3079 : // Displaced Constraints
3080 3667 : if (_fe_problem.getDisplacedProblem())
3081 1096 : constraintJacobians(jacobian_to_view, true);
3082 3667 : }
3083 : }
3084 527170 : 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 527170 : 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 527170 : if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
3096 517356 : 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 527170 : 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 430405 : 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 430405 : std::map<std::string, std::set<unsigned int>> bc_involved_vars;
3115 430405 : const std::set<BoundaryID> & all_boundary_ids = _mesh.getBoundaryIDs();
3116 2159970 : 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 1729565 : if (nbc_warehouse->hasActiveBoundaryObjects(bid))
3121 : {
3122 852003 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
3123 1789662 : for (const auto & bc : bcs)
3124 : {
3125 : const std::vector<MooseVariableFEBase *> & coupled_moose_vars =
3126 937659 : 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 937659 : std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3131 995931 : 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 937659 : 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 899375 : for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3144 468970 : _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 430405 : auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number());
3152 :
3153 : // Compute Jacobians for NodalBCBases
3154 430405 : const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange();
3155 22418237 : for (const auto & bnode : bnd_nodes)
3156 : {
3157 21987832 : BoundaryID boundary_id = bnode->_bnd_id;
3158 21987832 : Node * node = bnode->_node;
3159 :
3160 32263459 : if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
3161 10275627 : node->processor_id() == processor_id())
3162 : {
3163 8015271 : _fe_problem.reinitNodeFace(node, boundary_id, 0);
3164 :
3165 8015271 : const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
3166 17054291 : for (const auto & bc : bcs)
3167 : {
3168 : // Get the set of involved MOOSE vars for this BC
3169 9039020 : 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 22755976 : for (const auto & it : coupling_entries)
3176 : {
3177 13716956 : 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 13716956 : if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
3184 9128404 : bc->computeOffDiagJacobian(jvar);
3185 : }
3186 :
3187 9039020 : const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars();
3188 9039391 : 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 430405 : _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{});
3197 430405 : }
3198 : }
3199 527170 : PARALLEL_CATCH;
3200 :
3201 527170 : 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 527170 : if (_has_nodalbc_diag_save_in)
3206 23 : _fe_problem.getAuxiliarySystem().solution().close();
3207 :
3208 527170 : 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 527170 : _app.solutionInvalidity().syncIteration();
3214 527170 : _app.solutionInvalidity().solutionInvalidAccumulation();
3215 528586 : }
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 527890 : NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
3242 : {
3243 1583670 : TIME_SECTION("computeJacobianTags", 5);
3244 :
3245 527890 : FloatingPointExceptionGuard fpe_guard(_app);
3246 :
3247 : try
3248 : {
3249 527890 : 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 527882 : }
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 370685 : NonlinearSystemBase::updateActive(THREAD_ID tid)
3360 : {
3361 370685 : _element_dampers.updateActive(tid);
3362 370685 : _nodal_dampers.updateActive(tid);
3363 370685 : _integrated_bcs.updateActive(tid);
3364 370685 : _dg_kernels.updateActive(tid);
3365 370685 : _interface_kernels.updateActive(tid);
3366 370685 : _dirac_kernels.updateActive(tid);
3367 370685 : _kernels.updateActive(tid);
3368 370685 : _nodal_kernels.updateActive(tid);
3369 :
3370 370685 : if (tid == 0)
3371 : {
3372 340636 : _general_dampers.updateActive();
3373 340636 : _nodal_bcs.updateActive();
3374 340636 : _preset_nodal_bcs.updateActive();
3375 340636 : _ad_preset_nodal_bcs.updateActive();
3376 340636 : _splits.updateActive();
3377 340636 : _constraints.updateActive();
3378 340636 : _scalar_kernels.updateActive();
3379 :
3380 : #ifdef MOOSE_KOKKOS_ENABLED
3381 224572 : _kokkos_kernels.updateActive();
3382 224572 : _kokkos_nodal_kernels.updateActive();
3383 224572 : _kokkos_integrated_bcs.updateActive();
3384 224572 : _kokkos_nodal_bcs.updateActive();
3385 224572 : _kokkos_preset_nodal_bcs.updateActive();
3386 : #endif
3387 : }
3388 370685 : }
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 3866387 : NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
3478 : {
3479 3866387 : _fe_problem.clearDiracInfo();
3480 :
3481 3866387 : std::set<const Elem *> dirac_elements;
3482 :
3483 3866387 : if (_dirac_kernels.hasActiveObjects())
3484 : {
3485 192085 : TIME_SECTION("computeDirac", 3, "Computing DiracKernels");
3486 :
3487 : // TODO: Need a threading fix... but it's complicated!
3488 79882 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3489 : {
3490 41481 : const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
3491 96376 : for (const auto & dkernel : dkernels)
3492 : {
3493 54911 : dkernel->clearPoints();
3494 54911 : dkernel->addPoints();
3495 : }
3496 : }
3497 :
3498 38401 : ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
3499 :
3500 38401 : _fe_problem.getDiracElements(dirac_elements);
3501 :
3502 38401 : DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
3503 : // TODO: Make Dirac work thread!
3504 : // Threads::parallel_reduce(range, cd);
3505 :
3506 38401 : cd(range);
3507 38401 : }
3508 3866371 : }
3509 :
3510 : NumericVector<Number> &
3511 0 : NonlinearSystemBase::residualCopy()
3512 : {
3513 0 : if (!_residual_copy.get())
3514 0 : _residual_copy = NumericVector<Number>::build(_communicator);
3515 :
3516 0 : return *_residual_copy;
3517 : }
3518 :
3519 : NumericVector<Number> &
3520 360 : NonlinearSystemBase::residualGhosted()
3521 : {
3522 360 : _need_residual_ghosted = true;
3523 360 : if (!_residual_ghosted)
3524 : {
3525 : // The first time we realize we need a ghosted residual vector,
3526 : // we add it.
3527 576 : _residual_ghosted = &addVector("residual_ghosted", false, GHOSTED);
3528 :
3529 : // If we've already realized we need time and/or non-time
3530 : // residual vectors, but we haven't yet realized they need to be
3531 : // ghosted, fix that now.
3532 : //
3533 : // If an application changes its mind, the libMesh API lets us
3534 : // change the vector.
3535 288 : if (_Re_time)
3536 : {
3537 50 : const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
3538 50 : _Re_time = &system().add_vector(vector_name, false, GHOSTED);
3539 50 : }
3540 288 : if (_Re_non_time)
3541 : {
3542 288 : const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
3543 288 : _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
3544 288 : }
3545 : }
3546 360 : return *_residual_ghosted;
3547 : }
3548 :
3549 : void
3550 69961 : NonlinearSystemBase::augmentSparsity(SparsityPattern::Graph & sparsity,
3551 : std::vector<dof_id_type> & n_nz,
3552 : std::vector<dof_id_type> & n_oz)
3553 : {
3554 69961 : if (_add_implicit_geometric_coupling_entries_to_jacobian)
3555 : {
3556 1231 : _fe_problem.updateGeomSearch();
3557 :
3558 1231 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
3559 :
3560 1231 : findImplicitGeometricCouplingEntries(_fe_problem.geomSearchData(), graph);
3561 :
3562 1231 : if (_fe_problem.getDisplacedProblem())
3563 201 : findImplicitGeometricCouplingEntries(_fe_problem.getDisplacedProblem()->geomSearchData(),
3564 : graph);
3565 :
3566 1231 : const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
3567 1231 : const dof_id_type end_dof_on_proc = dofMap().end_dof(processor_id());
3568 :
3569 : // The total number of dofs on and off processor
3570 1231 : const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
3571 1231 : const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
3572 :
3573 2088 : for (const auto & git : graph)
3574 : {
3575 857 : dof_id_type dof = git.first;
3576 857 : dof_id_type local_dof = dof - first_dof_on_proc;
3577 :
3578 857 : if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
3579 126 : continue;
3580 :
3581 731 : const auto & row = git.second;
3582 :
3583 731 : SparsityPattern::Row & sparsity_row = sparsity[local_dof];
3584 :
3585 731 : unsigned int original_row_length = sparsity_row.size();
3586 :
3587 731 : sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
3588 :
3589 1462 : SparsityPattern::sort_row(
3590 731 : sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
3591 :
3592 : // Fix up nonzero arrays
3593 4807 : for (const auto & coupled_dof : row)
3594 : {
3595 4076 : if (coupled_dof < first_dof_on_proc || coupled_dof >= end_dof_on_proc)
3596 : {
3597 1224 : if (n_oz[local_dof] < n_dofs_not_on_proc)
3598 612 : n_oz[local_dof]++;
3599 : }
3600 : else
3601 : {
3602 3464 : if (n_nz[local_dof] < n_dofs_on_proc)
3603 3464 : n_nz[local_dof]++;
3604 : }
3605 : }
3606 : }
3607 1231 : }
3608 69961 : }
3609 :
3610 : void
3611 0 : NonlinearSystemBase::setSolutionUDot(const NumericVector<Number> & u_dot)
3612 : {
3613 0 : *_u_dot = u_dot;
3614 0 : }
3615 :
3616 : void
3617 0 : NonlinearSystemBase::setSolutionUDotDot(const NumericVector<Number> & u_dotdot)
3618 : {
3619 0 : *_u_dotdot = u_dotdot;
3620 0 : }
3621 :
3622 : void
3623 0 : NonlinearSystemBase::setSolutionUDotOld(const NumericVector<Number> & u_dot_old)
3624 : {
3625 0 : *_u_dot_old = u_dot_old;
3626 0 : }
3627 :
3628 : void
3629 0 : NonlinearSystemBase::setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old)
3630 : {
3631 0 : *_u_dotdot_old = u_dotdot_old;
3632 0 : }
3633 :
3634 : void
3635 15068 : NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
3636 : {
3637 15068 : if (_preconditioner.get() != nullptr)
3638 7 : mooseError("More than one active Preconditioner detected");
3639 :
3640 15061 : _preconditioner = pc;
3641 15061 : }
3642 :
3643 : MoosePreconditioner const *
3644 55666 : NonlinearSystemBase::getPreconditioner() const
3645 : {
3646 55666 : return _preconditioner.get();
3647 : }
3648 :
3649 : void
3650 322 : NonlinearSystemBase::setupDampers()
3651 : {
3652 322 : _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
3653 322 : }
3654 :
3655 : void
3656 1752 : NonlinearSystemBase::reinitIncrementAtQpsForDampers(THREAD_ID /*tid*/,
3657 : const std::set<MooseVariable *> & damped_vars)
3658 : {
3659 3526 : for (const auto & var : damped_vars)
3660 1774 : var->computeIncrementAtQps(*_increment_vec);
3661 1752 : }
3662 :
3663 : void
3664 7441 : NonlinearSystemBase::reinitIncrementAtNodeForDampers(THREAD_ID /*tid*/,
3665 : const std::set<MooseVariable *> & damped_vars)
3666 : {
3667 14882 : for (const auto & var : damped_vars)
3668 7441 : var->computeIncrementAtNode(*_increment_vec);
3669 7441 : }
3670 :
3671 : void
3672 43980 : NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
3673 : {
3674 : // Obtain all blocks and variables covered by all kernels
3675 43980 : std::set<SubdomainID> input_subdomains;
3676 43980 : std::set<std::string> kernel_variables;
3677 :
3678 43980 : bool global_kernels_exist = false;
3679 43980 : global_kernels_exist |= _scalar_kernels.hasActiveObjects();
3680 43980 : global_kernels_exist |= _nodal_kernels.hasActiveObjects();
3681 :
3682 43980 : _kernels.subdomainsCovered(input_subdomains, kernel_variables);
3683 43980 : _dg_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3684 43980 : _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3685 43980 : _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3686 43980 : _constraints.subdomainsCovered(input_subdomains, kernel_variables);
3687 :
3688 : #ifdef MOOSE_KOKKOS_ENABLED
3689 29469 : _kokkos_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3690 29469 : _kokkos_nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3691 : #endif
3692 :
3693 43980 : if (_fe_problem.haveFV())
3694 : {
3695 3538 : std::vector<FVElementalKernel *> fv_elemental_kernels;
3696 3538 : _fe_problem.theWarehouse()
3697 7076 : .query()
3698 3538 : .template condition<AttribSystem>("FVElementalKernel")
3699 3538 : .queryInto(fv_elemental_kernels);
3700 :
3701 7560 : for (auto fv_kernel : fv_elemental_kernels)
3702 : {
3703 4022 : if (fv_kernel->blockRestricted())
3704 1972 : for (auto block_id : fv_kernel->blockIDs())
3705 1103 : input_subdomains.insert(block_id);
3706 : else
3707 3153 : global_kernels_exist = true;
3708 4022 : kernel_variables.insert(fv_kernel->variable().name());
3709 :
3710 : // Check for lagrange multiplier
3711 4022 : if (dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel))
3712 268 : kernel_variables.insert(dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel)
3713 268 : ->lambdaVariable()
3714 134 : .name());
3715 : }
3716 :
3717 3538 : std::vector<FVFluxKernel *> fv_flux_kernels;
3718 3538 : _fe_problem.theWarehouse()
3719 7076 : .query()
3720 3538 : .template condition<AttribSystem>("FVFluxKernel")
3721 3538 : .queryInto(fv_flux_kernels);
3722 :
3723 8819 : for (auto fv_kernel : fv_flux_kernels)
3724 : {
3725 5281 : if (fv_kernel->blockRestricted())
3726 2350 : for (auto block_id : fv_kernel->blockIDs())
3727 1250 : input_subdomains.insert(block_id);
3728 : else
3729 4181 : global_kernels_exist = true;
3730 5281 : kernel_variables.insert(fv_kernel->variable().name());
3731 : }
3732 :
3733 3538 : std::vector<FVInterfaceKernel *> fv_interface_kernels;
3734 3538 : _fe_problem.theWarehouse()
3735 7076 : .query()
3736 3538 : .template condition<AttribSystem>("FVInterfaceKernel")
3737 3538 : .queryInto(fv_interface_kernels);
3738 :
3739 3843 : for (auto fvik : fv_interface_kernels)
3740 305 : if (auto scalar_fvik = dynamic_cast<FVScalarLagrangeMultiplierInterface *>(fvik))
3741 14 : kernel_variables.insert(scalar_fvik->lambdaVariable().name());
3742 :
3743 3538 : std::vector<FVFluxBC *> fv_flux_bcs;
3744 3538 : _fe_problem.theWarehouse()
3745 7076 : .query()
3746 3538 : .template condition<AttribSystem>("FVFluxBC")
3747 3538 : .queryInto(fv_flux_bcs);
3748 :
3749 5582 : for (auto fvbc : fv_flux_bcs)
3750 2044 : if (auto scalar_fvbc = dynamic_cast<FVBoundaryScalarLagrangeMultiplierConstraint *>(fvbc))
3751 14 : kernel_variables.insert(scalar_fvbc->lambdaVariable().name());
3752 3538 : }
3753 :
3754 : // Check kernel coverage of subdomains (blocks) in your mesh
3755 43980 : if (!global_kernels_exist)
3756 : {
3757 39941 : std::set<SubdomainID> difference;
3758 39941 : std::set_difference(mesh_subdomains.begin(),
3759 : mesh_subdomains.end(),
3760 : input_subdomains.begin(),
3761 : input_subdomains.end(),
3762 : std::inserter(difference, difference.end()));
3763 :
3764 : // there supposed to be no kernels on this lower-dimensional subdomain
3765 40179 : for (const auto & id : _mesh.interiorLowerDBlocks())
3766 238 : difference.erase(id);
3767 40179 : for (const auto & id : _mesh.boundaryLowerDBlocks())
3768 238 : difference.erase(id);
3769 :
3770 39941 : if (!difference.empty())
3771 : {
3772 : std::vector<SubdomainID> difference_vec =
3773 8 : std::vector<SubdomainID>(difference.begin(), difference.end());
3774 8 : std::vector<SubdomainName> difference_names = _mesh.getSubdomainNames(difference_vec);
3775 8 : std::stringstream missing_block_names;
3776 8 : std::copy(difference_names.begin(),
3777 : difference_names.end(),
3778 8 : std::ostream_iterator<std::string>(missing_block_names, " "));
3779 8 : std::stringstream missing_block_ids;
3780 8 : std::copy(difference.begin(),
3781 : difference.end(),
3782 8 : std::ostream_iterator<unsigned int>(missing_block_ids, " "));
3783 :
3784 8 : mooseError("Each subdomain must contain at least one Kernel.\nThe following block(s) lack an "
3785 8 : "active kernel: " +
3786 8 : missing_block_names.str(),
3787 : " (ids: ",
3788 8 : missing_block_ids.str(),
3789 : ")");
3790 0 : }
3791 39933 : }
3792 :
3793 : // Check kernel use of variables
3794 43972 : std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
3795 :
3796 43972 : std::set<VariableName> difference;
3797 43972 : std::set_difference(variables.begin(),
3798 : variables.end(),
3799 : kernel_variables.begin(),
3800 : kernel_variables.end(),
3801 : std::inserter(difference, difference.end()));
3802 :
3803 : // skip checks for varaibles defined on lower-dimensional subdomain
3804 43972 : std::set<VariableName> vars(difference);
3805 44335 : for (auto & var_name : vars)
3806 : {
3807 363 : auto blks = getSubdomainsForVar(var_name);
3808 744 : for (const auto & id : blks)
3809 381 : if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
3810 381 : difference.erase(var_name);
3811 363 : }
3812 :
3813 43972 : if (!difference.empty())
3814 : {
3815 8 : std::stringstream missing_kernel_vars;
3816 8 : std::copy(difference.begin(),
3817 : difference.end(),
3818 8 : std::ostream_iterator<std::string>(missing_kernel_vars, " "));
3819 8 : mooseError("Each variable must be referenced by at least one active Kernel.\nThe following "
3820 8 : "variable(s) lack an active kernel: " +
3821 8 : missing_kernel_vars.str());
3822 0 : }
3823 43964 : }
3824 :
3825 : bool
3826 29393 : NonlinearSystemBase::containsTimeKernel()
3827 : {
3828 29393 : auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3829 :
3830 29393 : return time_kernels.hasActiveObjects();
3831 : }
3832 :
3833 : std::vector<std::string>
3834 0 : NonlinearSystemBase::timeKernelVariableNames()
3835 : {
3836 0 : std::vector<std::string> variable_names;
3837 0 : const auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3838 0 : if (time_kernels.hasActiveObjects())
3839 0 : for (const auto & kernel : time_kernels.getObjects())
3840 0 : variable_names.push_back(kernel->variable().name());
3841 :
3842 0 : return variable_names;
3843 0 : }
3844 :
3845 : bool
3846 30339 : NonlinearSystemBase::needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3847 : {
3848 : // IntegratedBCs are for now the only objects we consider to be consuming
3849 : // matprops on boundaries.
3850 30339 : if (_integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid))
3851 6820 : for (const auto & bc : _integrated_bcs.getActiveBoundaryObjects(bnd_id, tid))
3852 4326 : if (std::static_pointer_cast<MaterialPropertyInterface>(bc)->getMaterialPropertyCalled())
3853 1604 : return true;
3854 :
3855 : // Thin layer heat transfer in the heat_transfer module is being used on a boundary even though
3856 : // it's an interface kernel. That boundary is external, on both sides of a gap in a mesh
3857 28735 : if (_interface_kernels.hasActiveBoundaryObjects(bnd_id, tid))
3858 391 : for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3859 311 : if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3860 231 : return true;
3861 :
3862 : // Because MortarConstraints do not inherit from BoundaryRestrictable, they are not sorted
3863 : // by boundary in the MooseObjectWarehouse. So for now, we return true for all boundaries
3864 : // Note: constraints are not threaded at this time
3865 28504 : if (_constraints.hasActiveObjects(/*tid*/ 0))
3866 2535 : for (const auto & ct : _constraints.getActiveObjects(/*tid*/ 0))
3867 2124 : if (auto mpi = std::dynamic_pointer_cast<MaterialPropertyInterface>(ct);
3868 2124 : mpi && mpi->getMaterialPropertyCalled())
3869 2124 : return true;
3870 27539 : return false;
3871 : }
3872 :
3873 : bool
3874 3238 : NonlinearSystemBase::needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
3875 : {
3876 : // InterfaceKernels are for now the only objects we consider to be consuming matprops on internal
3877 : // boundaries.
3878 3238 : if (_interface_kernels.hasActiveBoundaryObjects(bnd_id, tid))
3879 482 : for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3880 376 : if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3881 270 : return true;
3882 2968 : return false;
3883 : }
3884 :
3885 : bool
3886 13378 : NonlinearSystemBase::needInternalNeighborSideMaterial(SubdomainID subdomain_id, THREAD_ID tid) const
3887 : {
3888 : // DGKernels are for now the only objects we consider to be consuming matprops on
3889 : // internal sides.
3890 13378 : if (_dg_kernels.hasActiveBlockObjects(subdomain_id, tid))
3891 788 : for (const auto & dg : _dg_kernels.getActiveBlockObjects(subdomain_id, tid))
3892 654 : if (std::static_pointer_cast<MaterialPropertyInterface>(dg)->getMaterialPropertyCalled())
3893 495 : return true;
3894 : // NOTE:
3895 : // HDG kernels do not require face material properties on internal sides at this time.
3896 : // The idea is to have element locality of HDG for hybridization
3897 12883 : return false;
3898 : }
3899 :
3900 : bool
3901 0 : NonlinearSystemBase::doingDG() const
3902 : {
3903 0 : return _doing_dg;
3904 : }
3905 :
3906 : void
3907 556 : NonlinearSystemBase::setPreviousNewtonSolution(const NumericVector<Number> & soln)
3908 : {
3909 556 : if (hasVector(Moose::PREVIOUS_NL_SOLUTION_TAG))
3910 556 : getVector(Moose::PREVIOUS_NL_SOLUTION_TAG) = soln;
3911 556 : }
3912 :
3913 : void
3914 3870468 : NonlinearSystemBase::mortarConstraints(const Moose::ComputeType compute_type,
3915 : const std::set<TagID> & vector_tags,
3916 : const std::set<TagID> & matrix_tags)
3917 : {
3918 : parallel_object_only();
3919 :
3920 : try
3921 : {
3922 3879587 : for (auto & map_pr : _undisplaced_mortar_functors)
3923 9119 : map_pr.second(compute_type, vector_tags, matrix_tags);
3924 :
3925 3875156 : for (auto & map_pr : _displaced_mortar_functors)
3926 4688 : map_pr.second(compute_type, vector_tags, matrix_tags);
3927 : }
3928 0 : catch (MetaPhysicL::LogicError &)
3929 : {
3930 0 : mooseError(
3931 : "We caught a MetaPhysicL error in NonlinearSystemBase::mortarConstraints. This is very "
3932 : "likely due to AD not having a sufficiently large derivative container size. Please run "
3933 : "MOOSE configure with the '--with-derivative-size=<n>' option");
3934 0 : }
3935 3870468 : }
3936 :
3937 : void
3938 439 : NonlinearSystemBase::setupScalingData()
3939 : {
3940 439 : if (_auto_scaling_initd)
3941 0 : return;
3942 :
3943 : // Want the libMesh count of variables, not MOOSE, e.g. I don't care about array variable counts
3944 439 : const auto n_vars = system().n_vars();
3945 :
3946 439 : if (_scaling_group_variables.empty())
3947 : {
3948 427 : _var_to_group_var.reserve(n_vars);
3949 427 : _num_scaling_groups = n_vars;
3950 :
3951 1128 : for (const auto var_number : make_range(n_vars))
3952 701 : _var_to_group_var.emplace(var_number, var_number);
3953 : }
3954 : else
3955 : {
3956 12 : std::set<unsigned int> var_numbers, var_numbers_covered, var_numbers_not_covered;
3957 48 : for (const auto var_number : make_range(n_vars))
3958 36 : var_numbers.insert(var_number);
3959 :
3960 12 : _num_scaling_groups = _scaling_group_variables.size();
3961 :
3962 24 : for (const auto group_index : index_range(_scaling_group_variables))
3963 48 : for (const auto & var_name : _scaling_group_variables[group_index])
3964 : {
3965 36 : if (!hasVariable(var_name) && !hasScalarVariable(var_name))
3966 0 : mooseError("'",
3967 : var_name,
3968 : "', provided to the 'scaling_group_variables' parameter, does not exist in "
3969 : "the nonlinear system.");
3970 :
3971 : const MooseVariableBase & var =
3972 36 : hasVariable(var_name)
3973 24 : ? static_cast<MooseVariableBase &>(getVariable(0, var_name))
3974 60 : : static_cast<MooseVariableBase &>(getScalarVariable(0, var_name));
3975 36 : auto map_pair = _var_to_group_var.emplace(var.number(), group_index);
3976 36 : if (!map_pair.second)
3977 0 : mooseError("Variable ", var_name, " is contained in multiple scaling grouplings");
3978 36 : var_numbers_covered.insert(var.number());
3979 : }
3980 :
3981 12 : std::set_difference(var_numbers.begin(),
3982 : var_numbers.end(),
3983 : var_numbers_covered.begin(),
3984 : var_numbers_covered.end(),
3985 : std::inserter(var_numbers_not_covered, var_numbers_not_covered.begin()));
3986 :
3987 12 : _num_scaling_groups = _scaling_group_variables.size() + var_numbers_not_covered.size();
3988 :
3989 12 : auto index = static_cast<unsigned int>(_scaling_group_variables.size());
3990 12 : for (auto var_number : var_numbers_not_covered)
3991 0 : _var_to_group_var.emplace(var_number, index++);
3992 12 : }
3993 :
3994 439 : _variable_autoscaled.resize(n_vars, true);
3995 439 : const auto & number_to_var_map = _vars[0].numberToVariableMap();
3996 :
3997 439 : if (_ignore_variables_for_autoscaling.size())
3998 50 : for (const auto i : index_range(_variable_autoscaled))
3999 40 : if (std::find(_ignore_variables_for_autoscaling.begin(),
4000 : _ignore_variables_for_autoscaling.end(),
4001 80 : libmesh_map_find(number_to_var_map, i)->name()) !=
4002 80 : _ignore_variables_for_autoscaling.end())
4003 20 : _variable_autoscaled[i] = false;
4004 :
4005 439 : _auto_scaling_initd = true;
4006 : }
4007 :
4008 : bool
4009 1112 : NonlinearSystemBase::computeScaling()
4010 : {
4011 1112 : if (_compute_scaling_once && _computed_scaling)
4012 460 : return true;
4013 :
4014 652 : _console << "\nPerforming automatic scaling calculation\n" << std::endl;
4015 :
4016 3260 : TIME_SECTION("computeScaling", 3, "Computing Automatic Scaling");
4017 :
4018 : // It's funny but we need to assemble our vector of scaling factors here otherwise we will be
4019 : // applying scaling factors of 0 during Assembly of our scaling Jacobian
4020 652 : assembleScalingVector();
4021 :
4022 : // container for repeated access of element global dof indices
4023 652 : std::vector<dof_id_type> dof_indices;
4024 :
4025 652 : if (!_auto_scaling_initd)
4026 439 : setupScalingData();
4027 :
4028 652 : std::vector<Real> inverse_scaling_factors(_num_scaling_groups, 0);
4029 652 : std::vector<Real> resid_inverse_scaling_factors(_num_scaling_groups, 0);
4030 652 : std::vector<Real> jac_inverse_scaling_factors(_num_scaling_groups, 0);
4031 652 : auto & dof_map = dofMap();
4032 :
4033 : // what types of scaling do we want?
4034 652 : bool jac_scaling = _resid_vs_jac_scaling_param < 1. - TOLERANCE;
4035 652 : bool resid_scaling = _resid_vs_jac_scaling_param > TOLERANCE;
4036 :
4037 652 : const NumericVector<Number> & scaling_residual = RHS();
4038 :
4039 652 : if (jac_scaling)
4040 : {
4041 : // if (!_auto_scaling_initd)
4042 : // We need to reinit this when the number of dofs changes
4043 : // but there is no good way to track that
4044 : // In theory, it is the job of libmesh system to track this,
4045 : // but this special matrix is not owned by libMesh system
4046 : // Let us reinit eveytime since it is not expensive
4047 : {
4048 602 : auto init_vector = NumericVector<Number>::build(this->comm());
4049 602 : init_vector->init(system().n_dofs(), system().n_local_dofs(), /*fast=*/false, PARALLEL);
4050 :
4051 602 : _scaling_matrix->clear();
4052 602 : _scaling_matrix->init(*init_vector);
4053 602 : }
4054 :
4055 602 : _fe_problem.computingScalingJacobian(true);
4056 : // Dispatch to derived classes to ensure that we use the correct matrix tag
4057 602 : computeScalingJacobian();
4058 602 : _fe_problem.computingScalingJacobian(false);
4059 : }
4060 :
4061 652 : if (resid_scaling)
4062 : {
4063 50 : _fe_problem.computingScalingResidual(true);
4064 50 : _fe_problem.computingNonlinearResid(true);
4065 : // Dispatch to derived classes to ensure that we use the correct vector tag
4066 50 : computeScalingResidual();
4067 50 : _fe_problem.computingNonlinearResid(false);
4068 50 : _fe_problem.computingScalingResidual(false);
4069 : }
4070 :
4071 : // Did something bad happen during residual/Jacobian scaling computation?
4072 652 : if (_fe_problem.getFailNextNonlinearConvergenceCheck())
4073 0 : return false;
4074 :
4075 228330 : auto examine_dof_indices = [this,
4076 : jac_scaling,
4077 : resid_scaling,
4078 : &dof_map,
4079 : &jac_inverse_scaling_factors,
4080 : &resid_inverse_scaling_factors,
4081 : &scaling_residual](const auto & dof_indices, const auto var_number)
4082 : {
4083 1136309 : for (auto dof_index : dof_indices)
4084 907979 : if (dof_map.local_index(dof_index))
4085 : {
4086 903220 : if (jac_scaling)
4087 : {
4088 : // For now we will use the diagonal for determining scaling
4089 899103 : auto mat_value = (*_scaling_matrix)(dof_index, dof_index);
4090 899103 : auto & factor = jac_inverse_scaling_factors[_var_to_group_var[var_number]];
4091 899103 : factor = std::max(factor, std::abs(mat_value));
4092 : }
4093 903220 : if (resid_scaling)
4094 : {
4095 4117 : auto vec_value = scaling_residual(dof_index);
4096 4117 : auto & factor = resid_inverse_scaling_factors[_var_to_group_var[var_number]];
4097 4117 : factor = std::max(factor, std::abs(vec_value));
4098 : }
4099 : }
4100 228330 : };
4101 :
4102 : // Compute our scaling factors for the spatial field variables
4103 82245 : for (const auto & elem : _fe_problem.getCurrentAlgebraicElementRange())
4104 310836 : for (const auto i : make_range(system().n_vars()))
4105 229243 : if (_variable_autoscaled[i] && system().variable_type(i).family != SCALAR)
4106 : {
4107 228237 : dof_map.dof_indices(elem, dof_indices, i);
4108 228237 : examine_dof_indices(dof_indices, i);
4109 : }
4110 :
4111 1712 : for (const auto i : make_range(system().n_vars()))
4112 1060 : if (_variable_autoscaled[i] && system().variable_type(i).family == SCALAR)
4113 : {
4114 93 : dof_map.SCALAR_dof_indices(dof_indices, i);
4115 93 : examine_dof_indices(dof_indices, i);
4116 : }
4117 :
4118 652 : if (resid_scaling)
4119 50 : _communicator.max(resid_inverse_scaling_factors);
4120 652 : if (jac_scaling)
4121 602 : _communicator.max(jac_inverse_scaling_factors);
4122 :
4123 652 : if (jac_scaling && resid_scaling)
4124 0 : for (MooseIndex(inverse_scaling_factors) i = 0; i < inverse_scaling_factors.size(); ++i)
4125 : {
4126 : // Be careful not to take log(0)
4127 0 : if (!resid_inverse_scaling_factors[i])
4128 : {
4129 0 : if (!jac_inverse_scaling_factors[i])
4130 0 : inverse_scaling_factors[i] = 1;
4131 : else
4132 0 : inverse_scaling_factors[i] = jac_inverse_scaling_factors[i];
4133 : }
4134 0 : else if (!jac_inverse_scaling_factors[i])
4135 : // We know the resid is not zero
4136 0 : inverse_scaling_factors[i] = resid_inverse_scaling_factors[i];
4137 : else
4138 0 : inverse_scaling_factors[i] =
4139 0 : std::exp(_resid_vs_jac_scaling_param * std::log(resid_inverse_scaling_factors[i]) +
4140 0 : (1 - _resid_vs_jac_scaling_param) * std::log(jac_inverse_scaling_factors[i]));
4141 0 : }
4142 652 : else if (jac_scaling)
4143 602 : inverse_scaling_factors = jac_inverse_scaling_factors;
4144 50 : else if (resid_scaling)
4145 50 : inverse_scaling_factors = resid_inverse_scaling_factors;
4146 : else
4147 0 : mooseError("We shouldn't be calling this routine if we're not performing any scaling");
4148 :
4149 : // We have to make sure that our scaling values are not zero
4150 1688 : for (auto & scaling_factor : inverse_scaling_factors)
4151 1036 : if (scaling_factor == 0)
4152 40 : scaling_factor = 1;
4153 :
4154 : // Now flatten the group scaling factors to the individual variable scaling factors
4155 652 : std::vector<Real> flattened_inverse_scaling_factors(system().n_vars());
4156 1712 : for (const auto i : index_range(flattened_inverse_scaling_factors))
4157 1060 : flattened_inverse_scaling_factors[i] = inverse_scaling_factors[_var_to_group_var[i]];
4158 :
4159 : // Now set the scaling factors for the variables
4160 652 : applyScalingFactors(flattened_inverse_scaling_factors);
4161 652 : if (auto displaced_problem = _fe_problem.getDisplacedProblem().get())
4162 63 : displaced_problem->systemBaseNonlinear(number()).applyScalingFactors(
4163 : flattened_inverse_scaling_factors);
4164 :
4165 652 : _computed_scaling = true;
4166 652 : return true;
4167 652 : }
4168 :
4169 : void
4170 318548 : NonlinearSystemBase::assembleScalingVector()
4171 : {
4172 955644 : if (!hasVector("scaling_factors"))
4173 : // No variables have indicated they need scaling
4174 317998 : return;
4175 :
4176 1100 : auto & scaling_vector = getVector("scaling_factors");
4177 :
4178 550 : const auto & lm_mesh = _mesh.getMesh();
4179 550 : const auto & dof_map = dofMap();
4180 :
4181 550 : const auto & field_variables = _vars[0].fieldVariables();
4182 550 : const auto & scalar_variables = _vars[0].scalars();
4183 :
4184 550 : std::vector<dof_id_type> dof_indices;
4185 :
4186 550 : for (const Elem * const elem :
4187 13672 : as_range(lm_mesh.active_local_elements_begin(), lm_mesh.active_local_elements_end()))
4188 31020 : for (const auto * const field_var : field_variables)
4189 : {
4190 18448 : const auto & factors = field_var->arrayScalingFactor();
4191 36896 : for (const auto i : make_range(field_var->count()))
4192 : {
4193 18448 : dof_map.dof_indices(elem, dof_indices, field_var->number() + i);
4194 76448 : for (const auto dof : dof_indices)
4195 58000 : scaling_vector.set(dof, factors[i]);
4196 : }
4197 550 : }
4198 :
4199 662 : for (const auto * const scalar_var : scalar_variables)
4200 : {
4201 : mooseAssert(scalar_var->count() == 1,
4202 : "Scalar variables should always have only one component.");
4203 112 : dof_map.SCALAR_dof_indices(dof_indices, scalar_var->number());
4204 224 : for (const auto dof : dof_indices)
4205 112 : scaling_vector.set(dof, scalar_var->scalingFactor());
4206 : }
4207 :
4208 : // Parallel assemble
4209 550 : scaling_vector.close();
4210 :
4211 550 : if (auto * displaced_problem = _fe_problem.getDisplacedProblem().get())
4212 : // copy into the corresponding displaced system vector because they should be the exact same
4213 0 : displaced_problem->systemBaseNonlinear(number()).getVector("scaling_factors") = scaling_vector;
4214 550 : }
4215 :
4216 : bool
4217 317896 : NonlinearSystemBase::preSolve()
4218 : {
4219 : // Clear the iteration counters
4220 317896 : _current_l_its.clear();
4221 317896 : _current_nl_its = 0;
4222 :
4223 : // Initialize the solution vector using a predictor and known values from nodal bcs
4224 317896 : setInitialSolution();
4225 :
4226 : // Now that the initial solution has ben set, potentially perform a residual/Jacobian evaluation
4227 : // to determine variable scaling factors
4228 317896 : if (_automatic_scaling)
4229 : {
4230 1112 : const bool scaling_succeeded = computeScaling();
4231 1112 : if (!scaling_succeeded)
4232 0 : return false;
4233 : }
4234 :
4235 : // We do not know a priori what variable a global degree of freedom corresponds to, so we need a
4236 : // map from global dof to scaling factor. We just use a ghosted NumericVector for that mapping
4237 317896 : assembleScalingVector();
4238 :
4239 317896 : return true;
4240 : }
4241 :
4242 : void
4243 5612 : NonlinearSystemBase::destroyColoring()
4244 : {
4245 5612 : if (matrixFromColoring())
4246 117 : LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
4247 5612 : }
4248 :
4249 : FieldSplitPreconditionerBase &
4250 0 : NonlinearSystemBase::getFieldSplitPreconditioner()
4251 : {
4252 0 : if (!_fsp)
4253 0 : mooseError("No field split preconditioner is present for this system");
4254 :
4255 0 : return *_fsp;
4256 : }
|
