https://mooseframework.inl.gov
NonlinearSystemBase.C
Go to the documentation of this file.
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"
23 #include "ComputeFVFluxThread.h"
24 #include "ComputeJacobianThread.h"
28 #include "ComputeDiracThread.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"
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"
69 #include "ADKernel.h"
70 #include "ADDirichletBCBase.h"
71 #include "Moose.h"
72 #include "ConsoleStream.h"
73 #include "MooseError.h"
74 #include "FVElementalKernel.h"
77 #include "FVFluxKernel.h"
79 #include "UserObject.h"
81 #include "HDGKernel.h"
82 
83 // libMesh
84 #include "libmesh/nonlinear_solver.h"
85 #include "libmesh/quadrature_gauss.h"
86 #include "libmesh/dense_vector.h"
87 #include "libmesh/boundary_info.h"
88 #include "libmesh/petsc_matrix.h"
89 #include "libmesh/petsc_vector.h"
90 #include "libmesh/petsc_nonlinear_solver.h"
91 #include "libmesh/numeric_vector.h"
92 #include "libmesh/mesh.h"
93 #include "libmesh/dense_subvector.h"
94 #include "libmesh/dense_submatrix.h"
95 #include "libmesh/dof_map.h"
96 #include "libmesh/sparse_matrix.h"
97 #include "libmesh/petsc_matrix.h"
98 #include "libmesh/default_coupling.h"
99 #include "libmesh/diagonal_matrix.h"
100 #include "libmesh/fe_interface.h"
101 #include "libmesh/petsc_solver_exception.h"
102 
103 #include <ios>
104 
105 #include "petscsnes.h"
106 #include <PetscDMMoose.h>
107 EXTERN_C_BEGIN
108 extern PetscErrorCode DMCreate_Moose(DM);
109 EXTERN_C_END
110 
111 using namespace libMesh;
112 
114  System & sys,
115  const std::string & name)
116  : SolverSystem(fe_problem, fe_problem, name, Moose::VAR_SOLVER),
117  PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
118  _sys(sys),
119  _last_nl_rnorm(0.),
120  _current_nl_its(0),
121  _residual_ghosted(NULL),
122  _Re_time_tag(-1),
123  _Re_time(NULL),
124  _Re_non_time_tag(-1),
125  _Re_non_time(NULL),
126  _scalar_kernels(/*threaded=*/false),
127  _nodal_bcs(/*threaded=*/false),
128  _preset_nodal_bcs(/*threaded=*/false),
129  _ad_preset_nodal_bcs(/*threaded=*/false),
130 #ifdef MOOSE_KOKKOS_ENABLED
131  _kokkos_kernels(/*threaded=*/false),
132  _kokkos_integrated_bcs(/*threaded=*/false),
133  _kokkos_nodal_bcs(/*threaded=*/false),
134  _kokkos_preset_nodal_bcs(/*threaded=*/false),
135  _kokkos_nodal_kernels(/*threaded=*/false),
136 #endif
137  _general_dampers(/*threaded=*/false),
138  _splits(/*threaded=*/false),
139  _increment_vec(NULL),
140  _use_finite_differenced_preconditioner(false),
141  _fdcoloring(nullptr),
142  _fsp(nullptr),
143  _add_implicit_geometric_coupling_entries_to_jacobian(false),
144  _assemble_constraints_separately(false),
145  _need_residual_ghosted(false),
146  _debugging_residuals(false),
147  _doing_dg(false),
148  _n_iters(0),
149  _n_linear_iters(0),
150  _n_residual_evaluations(0),
151  _final_residual(0.),
152  _computing_pre_smo_residual(false),
153  _pre_smo_residual(0),
154  _initial_residual(0),
155  _use_pre_smo_residual(false),
156  _print_all_var_norms(false),
157  _has_save_in(false),
158  _has_diag_save_in(false),
159  _has_nodalbc_save_in(false),
160  _has_nodalbc_diag_save_in(false),
161  _computed_scaling(false),
162  _compute_scaling_once(true),
163  _resid_vs_jac_scaling_param(0),
164  _off_diagonals_in_auto_scaling(false),
165  _auto_scaling_initd(false)
166 {
168  // Don't need to add the matrix - it already exists (for now)
170 
171  // The time matrix tag is not normally used - but must be added to the system
172  // in case it is so that objects can have 'time' in their matrix tags by default
173  _fe_problem.addMatrixTag("TIME");
174 
175  _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
176 
178 
180  {
181  auto & dof_map = _sys.get_dof_map();
182  dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting());
183  dof_map.set_implicit_neighbor_dofs(false);
184  }
185 }
186 
188 
189 void
191 {
193 
194  if (_fe_problem.hasDampers())
195  setupDampers();
196 
197  if (_residual_copy.get())
198  _residual_copy->init(_sys.n_dofs(), false, SERIAL);
199 
200 #ifdef MOOSE_KOKKOS_ENABLED
203 #endif
204 }
205 
206 void
208 {
209  // reinit is called on meshChanged() in FEProblemBase. We could implement meshChanged() instead.
210  // Subdomains might have changed
211  for (auto & functor : _displaced_mortar_functors)
212  functor.second.setupMortarMaterials();
213  for (auto & functor : _undisplaced_mortar_functors)
214  functor.second.setupMortarMaterials();
215 }
216 
217 void
219 {
221  nonlinearSolver()->jacobian = NULL;
222 }
223 
224 void
226 {
227  TIME_SECTION("nlInitialSetup", 2, "Setting Up Nonlinear System");
228 
230 
231  {
232  TIME_SECTION("kernelsInitialSetup", 2, "Setting Up Kernels/BCs/Constraints");
233 
234  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
235  {
236  _kernels.initialSetup(tid);
239  if (_doing_dg)
242 
246 
247  if (_fe_problem.haveFV())
248  {
249  std::vector<FVElementalKernel *> fv_elemental_kernels;
251  .query()
252  .template condition<AttribSystem>("FVElementalKernel")
253  .template condition<AttribThread>(tid)
254  .queryInto(fv_elemental_kernels);
255 
256  for (auto * fv_kernel : fv_elemental_kernels)
257  fv_kernel->initialSetup();
258 
259  std::vector<FVFluxKernel *> fv_flux_kernels;
261  .query()
262  .template condition<AttribSystem>("FVFluxKernel")
263  .template condition<AttribThread>(tid)
264  .queryInto(fv_flux_kernels);
265 
266  for (auto * fv_kernel : fv_flux_kernels)
267  fv_kernel->initialSetup();
268  }
269  }
270 
277 
278 #ifdef MOOSE_KOKKOS_ENABLED
283 #endif
284  }
285 
286  {
287  TIME_SECTION("mortarSetup", 2, "Initializing Mortar Interfaces");
288 
289  auto create_mortar_functors = [this](const bool displaced)
290  {
291  // go over mortar interfaces and construct functors
292  const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
293  for (const auto & mortar_interface : mortar_interfaces)
294  {
295  const auto primary_secondary_boundary_pair = mortar_interface.first;
296  if (!_constraints.hasActiveMortarConstraints(primary_secondary_boundary_pair, displaced))
297  continue;
298 
299  const auto & mortar_generation_object = mortar_interface.second;
300 
301  auto & mortar_constraints =
302  _constraints.getActiveMortarConstraints(primary_secondary_boundary_pair, displaced);
303 
304  auto & subproblem = displaced
305  ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
306  : static_cast<SubProblem &>(_fe_problem);
307 
308  auto & mortar_functors =
310 
311  mortar_functors.emplace(primary_secondary_boundary_pair,
312  ComputeMortarFunctor(mortar_constraints,
313  mortar_generation_object,
314  subproblem,
315  _fe_problem,
316  displaced,
317  subproblem.assembly(0, number())));
318  }
319  };
320 
321  create_mortar_functors(false);
322  create_mortar_functors(true);
323  }
324 
325  if (_automatic_scaling)
326  {
328  _scaling_matrix = std::make_unique<OffDiagonalScalingMatrix<Number>>(_communicator);
329  else
330  _scaling_matrix = std::make_unique<DiagonalMatrix<Number>>(_communicator);
331  }
332 
333  if (_preconditioner)
334  _preconditioner->initialSetup();
335 }
336 
337 void
339 {
341 
342  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
343  {
344  _kernels.timestepSetup(tid);
347  if (_doing_dg)
353 
354  if (_fe_problem.haveFV())
355  {
356  std::vector<FVFluxBC *> bcs;
358  .query()
359  .template condition<AttribSystem>("FVFluxBC")
360  .template condition<AttribThread>(tid)
361  .queryInto(bcs);
362 
363  std::vector<FVInterfaceKernel *> iks;
365  .query()
366  .template condition<AttribSystem>("FVInterfaceKernel")
367  .template condition<AttribThread>(tid)
368  .queryInto(iks);
369 
370  std::vector<FVFluxKernel *> kernels;
372  .query()
373  .template condition<AttribSystem>("FVFluxKernel")
374  .template condition<AttribThread>(tid)
375  .queryInto(kernels);
376 
377  for (auto * bc : bcs)
378  bc->timestepSetup();
379  for (auto * ik : iks)
380  ik->timestepSetup();
381  for (auto * kernel : kernels)
382  kernel->timestepSetup();
383  }
384  }
391 
392 #ifdef MOOSE_KOKKOS_ENABLED
397 #endif
398 }
399 
400 void
402 {
403  SolverSystem::customSetup(exec_type);
404 
405  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
406  {
407  _kernels.customSetup(exec_type, tid);
408  _nodal_kernels.customSetup(exec_type, tid);
409  _dirac_kernels.customSetup(exec_type, tid);
410  if (_doing_dg)
411  _dg_kernels.customSetup(exec_type, tid);
412  _interface_kernels.customSetup(exec_type, tid);
413  _element_dampers.customSetup(exec_type, tid);
414  _nodal_dampers.customSetup(exec_type, tid);
415  _integrated_bcs.customSetup(exec_type, tid);
416 
417  if (_fe_problem.haveFV())
418  {
419  std::vector<FVFluxBC *> bcs;
421  .query()
422  .template condition<AttribSystem>("FVFluxBC")
423  .template condition<AttribThread>(tid)
424  .queryInto(bcs);
425 
426  std::vector<FVInterfaceKernel *> iks;
428  .query()
429  .template condition<AttribSystem>("FVInterfaceKernel")
430  .template condition<AttribThread>(tid)
431  .queryInto(iks);
432 
433  std::vector<FVFluxKernel *> kernels;
435  .query()
436  .template condition<AttribSystem>("FVFluxKernel")
437  .template condition<AttribThread>(tid)
438  .queryInto(kernels);
439 
440  for (auto * bc : bcs)
441  bc->customSetup(exec_type);
442  for (auto * ik : iks)
443  ik->customSetup(exec_type);
444  for (auto * kernel : kernels)
445  kernel->customSetup(exec_type);
446  }
447  }
448  _scalar_kernels.customSetup(exec_type);
449  _constraints.customSetup(exec_type);
450  _general_dampers.customSetup(exec_type);
451  _nodal_bcs.customSetup(exec_type);
452  _preset_nodal_bcs.customSetup(exec_type);
454 
455 #ifdef MOOSE_KOKKOS_ENABLED
456  _kokkos_kernels.customSetup(exec_type);
459  _kokkos_nodal_bcs.customSetup(exec_type);
460 #endif
461 }
462 
463 void
465 {
466  if (_fsp)
467  _fsp->setupDM();
468 }
469 
470 void
471 NonlinearSystemBase::addKernel(const std::string & kernel_name,
472  const std::string & name,
473  InputParameters & parameters)
474 {
475  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
476  {
477  // Create the kernel object via the factory and add to warehouse
478  std::shared_ptr<KernelBase> kernel =
479  _factory.create<KernelBase>(kernel_name, name, parameters, tid);
480  _kernels.addObject(kernel, tid);
481  postAddResidualObject(*kernel);
482  // Add to theWarehouse, a centralized storage for all moose objects
483  _fe_problem.theWarehouse().add(kernel);
484  }
485 
486  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
487  _has_save_in = true;
488  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
489  _has_diag_save_in = true;
490 }
491 
492 void
493 NonlinearSystemBase::addHDGKernel(const std::string & kernel_name,
494  const std::string & name,
495  InputParameters & parameters)
496 {
497  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
498  {
499  // Create the kernel object via the factory and add to warehouse
500  auto kernel = _factory.create<HDGKernel>(kernel_name, name, parameters, tid);
501  _kernels.addObject(kernel, tid);
502  _hybridized_kernels.addObject(kernel, tid);
503  // Add to theWarehouse, a centralized storage for all moose objects
504  _fe_problem.theWarehouse().add(kernel);
505  postAddResidualObject(*kernel);
506  }
507 }
508 
509 void
510 NonlinearSystemBase::addNodalKernel(const std::string & kernel_name,
511  const std::string & name,
512  InputParameters & parameters)
513 {
514  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
515  {
516  // Create the kernel object via the factory and add to the warehouse
517  std::shared_ptr<NodalKernelBase> kernel =
518  _factory.create<NodalKernelBase>(kernel_name, name, parameters, tid);
519  _nodal_kernels.addObject(kernel, tid);
520  // Add to theWarehouse, a centralized storage for all moose objects
521  _fe_problem.theWarehouse().add(kernel);
522  postAddResidualObject(*kernel);
523  }
524 
525  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
526  _has_save_in = true;
527  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
528  _has_diag_save_in = true;
529 }
530 
531 void
532 NonlinearSystemBase::addScalarKernel(const std::string & kernel_name,
533  const std::string & name,
534  InputParameters & parameters)
535 {
536  std::shared_ptr<ScalarKernelBase> kernel =
537  _factory.create<ScalarKernelBase>(kernel_name, name, parameters);
538  postAddResidualObject(*kernel);
539  // Add to theWarehouse, a centralized storage for all moose objects
540  _fe_problem.theWarehouse().add(kernel);
541  _scalar_kernels.addObject(kernel);
542 }
543 
544 void
545 NonlinearSystemBase::addBoundaryCondition(const std::string & bc_name,
546  const std::string & name,
547  InputParameters & parameters)
548 {
549  // ThreadID
550  THREAD_ID tid = 0;
551 
552  // Create the object
553  std::shared_ptr<BoundaryCondition> bc =
554  _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
556 
557  // Active BoundaryIDs for the object
558  const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
559  auto bc_var = dynamic_cast<const MooseVariableFieldBase *>(&bc->variable());
560  _vars[tid].addBoundaryVar(boundary_ids, bc_var);
561 
562  // Cast to the various types of BCs
563  std::shared_ptr<NodalBCBase> nbc = std::dynamic_pointer_cast<NodalBCBase>(bc);
564  std::shared_ptr<IntegratedBCBase> ibc = std::dynamic_pointer_cast<IntegratedBCBase>(bc);
565 
566  // NodalBCBase
567  if (nbc)
568  {
569  if (nbc->checkNodalVar() && !nbc->variable().isNodal())
570  mooseError("Trying to use nodal boundary condition '",
571  nbc->name(),
572  "' on a non-nodal variable '",
573  nbc->variable().name(),
574  "'.");
575 
576  _nodal_bcs.addObject(nbc);
577  // Add to theWarehouse, a centralized storage for all moose objects
579  _vars[tid].addBoundaryVars(boundary_ids, nbc->getCoupledVars());
580 
581  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
582  _has_nodalbc_save_in = true;
583  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
585 
586  // DirichletBCs that are preset
587  std::shared_ptr<DirichletBCBase> dbc = std::dynamic_pointer_cast<DirichletBCBase>(bc);
588  if (dbc && dbc->preset())
590 
591  std::shared_ptr<ADDirichletBCBase> addbc = std::dynamic_pointer_cast<ADDirichletBCBase>(bc);
592  if (addbc && addbc->preset())
594  }
595 
596  // IntegratedBCBase
597  else if (ibc)
598  {
599  _integrated_bcs.addObject(ibc, tid);
600  // Add to theWarehouse, a centralized storage for all moose objects
602  _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
603 
604  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
605  _has_save_in = true;
606  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
607  _has_diag_save_in = true;
608 
609  for (tid = 1; tid < libMesh::n_threads(); tid++)
610  {
611  // Create the object
612  bc = _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
613 
614  // Give users opportunity to set some parameters
616 
617  // Active BoundaryIDs for the object
618  const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
619  _vars[tid].addBoundaryVar(boundary_ids, bc_var);
620 
621  ibc = std::static_pointer_cast<IntegratedBCBase>(bc);
622 
623  _integrated_bcs.addObject(ibc, tid);
624  _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
625  }
626  }
627 
628  else
629  mooseError("Unknown BoundaryCondition type for object named ", bc->name());
630 }
631 
632 void
633 NonlinearSystemBase::addConstraint(const std::string & c_name,
634  const std::string & name,
635  InputParameters & parameters)
636 {
637  std::shared_ptr<Constraint> constraint = _factory.create<Constraint>(c_name, name, parameters);
638  _constraints.addObject(constraint);
639  postAddResidualObject(*constraint);
640 
641  if (constraint && constraint->addCouplingEntriesToJacobian())
643 }
644 
645 void
646 NonlinearSystemBase::addDiracKernel(const std::string & kernel_name,
647  const std::string & name,
648  InputParameters & parameters)
649 {
650  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
651  {
652  std::shared_ptr<DiracKernelBase> kernel =
653  _factory.create<DiracKernelBase>(kernel_name, name, parameters, tid);
654  postAddResidualObject(*kernel);
655  _dirac_kernels.addObject(kernel, tid);
656  // Add to theWarehouse, a centralized storage for all moose objects
657  _fe_problem.theWarehouse().add(kernel);
658  }
659 }
660 
661 void
662 NonlinearSystemBase::addDGKernel(std::string dg_kernel_name,
663  const std::string & name,
664  InputParameters & parameters)
665 {
666  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
667  {
668  auto dg_kernel = _factory.create<DGKernelBase>(dg_kernel_name, name, parameters, tid);
669  _dg_kernels.addObject(dg_kernel, tid);
670  // Add to theWarehouse, a centralized storage for all moose objects
671  _fe_problem.theWarehouse().add(dg_kernel);
672  postAddResidualObject(*dg_kernel);
673  }
674 
675  _doing_dg = true;
676 
677  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
678  _has_save_in = true;
679  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
680  _has_diag_save_in = true;
681 }
682 
683 void
684 NonlinearSystemBase::addInterfaceKernel(std::string interface_kernel_name,
685  const std::string & name,
686  InputParameters & parameters)
687 {
688  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
689  {
690  std::shared_ptr<InterfaceKernelBase> interface_kernel =
691  _factory.create<InterfaceKernelBase>(interface_kernel_name, name, parameters, tid);
692  postAddResidualObject(*interface_kernel);
693 
694  const std::set<BoundaryID> & boundary_ids = interface_kernel->boundaryIDs();
695  auto ik_var = dynamic_cast<const MooseVariableFieldBase *>(&interface_kernel->variable());
696  _vars[tid].addBoundaryVar(boundary_ids, ik_var);
697 
698  _interface_kernels.addObject(interface_kernel, tid);
699  // Add to theWarehouse, a centralized storage for all moose objects
700  _fe_problem.theWarehouse().add(interface_kernel);
701  _vars[tid].addBoundaryVars(boundary_ids, interface_kernel->getCoupledVars());
702  }
703 }
704 
705 void
706 NonlinearSystemBase::addDamper(const std::string & damper_name,
707  const std::string & name,
708  InputParameters & parameters)
709 {
710  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
711  {
712  std::shared_ptr<Damper> damper = _factory.create<Damper>(damper_name, name, parameters, tid);
713 
714  // Attempt to cast to the damper types
715  std::shared_ptr<ElementDamper> ed = std::dynamic_pointer_cast<ElementDamper>(damper);
716  std::shared_ptr<NodalDamper> nd = std::dynamic_pointer_cast<NodalDamper>(damper);
717  std::shared_ptr<GeneralDamper> gd = std::dynamic_pointer_cast<GeneralDamper>(damper);
718 
719  if (gd)
720  {
722  break; // not threaded
723  }
724  else if (ed)
725  _element_dampers.addObject(ed, tid);
726  else if (nd)
727  _nodal_dampers.addObject(nd, tid);
728  else
729  mooseError("Invalid damper type");
730  }
731 }
732 
733 void
734 NonlinearSystemBase::addSplit(const std::string & split_name,
735  const std::string & name,
736  InputParameters & parameters)
737 {
738  std::shared_ptr<Split> split = _factory.create<Split>(split_name, name, parameters);
740  // Add to theWarehouse, a centralized storage for all moose objects
742 }
743 
744 std::shared_ptr<Split>
745 NonlinearSystemBase::getSplit(const std::string & name)
746 {
747  return _splits.getActiveObject(name);
748 }
749 
750 bool
752 {
754  return false;
755 
756  // The legacy behavior (#10464) _always_ performs the pre-SMO residual evaluation
757  // regardless of whether it is needed.
758  //
759  // This is not ideal and has been fixed by #23472. This legacy option ensures a smooth transition
760  // to the new behavior. Modules and Apps that want to migrate to the new behavior should set this
761  // parameter to false.
762  if (_app.parameters().get<bool>("use_legacy_initial_residual_evaluation_behavior"))
763  return true;
764 
765  return _use_pre_smo_residual;
766 }
767 
768 Real
770 {
772 }
773 
774 Real
776 {
778  mooseError("pre-SMO residual is requested but not evaluated.");
779 
780  return _pre_smo_residual;
781 }
782 
783 Real
785 {
786  return _initial_residual;
787 }
788 
789 void
791 {
792  _initial_residual = r;
793 }
794 
795 void
796 NonlinearSystemBase::zeroVectorForResidual(const std::string & vector_name)
797 {
798  for (unsigned int i = 0; i < _vecs_to_zero_for_residual.size(); ++i)
799  if (vector_name == _vecs_to_zero_for_residual[i])
800  return;
801 
802  _vecs_to_zero_for_residual.push_back(vector_name);
803 }
804 
805 void
807 {
808  _nl_vector_tags.clear();
809  _nl_vector_tags.insert(tag_id);
811 
813 
815 
817 }
818 
819 void
821 {
822  mooseDeprecated(" Please use computeResidualTag");
823 
824  computeResidualTag(residual, tag_id);
825 }
826 
827 void
828 NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
829 {
830  parallel_object_only();
831 
832  TIME_SECTION("nl::computeResidualTags", 5);
833 
836 
837  bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
838 
840 
841  // not suppose to do anythin on matrix
843 
845 
846  for (const auto & numeric_vec : _vecs_to_zero_for_residual)
847  if (hasVector(numeric_vec))
848  {
849  NumericVector<Number> & vec = getVector(numeric_vec);
850  vec.close();
851  vec.zero();
852  }
853 
854  try
855  {
856  zeroTaggedVectors(tags);
858  closeTaggedVectors(tags);
859 
860  if (required_residual)
861  {
862  auto & residual = getVector(residualVectorTag());
863  if (!_time_integrators.empty())
864  {
865  for (auto & ti : _time_integrators)
866  ti->postResidual(residual);
867  }
868  else
869  residual += *_Re_non_time;
870  residual.close();
871  }
873  // We don't want to do nodal bcs or anything else
874  return;
875 
876  computeNodalBCs(tags);
877  closeTaggedVectors(tags);
878 
879  // If we are debugging residuals we need one more assignment to have the ghosted copy up to
880  // date
881  if (_need_residual_ghosted && _debugging_residuals && required_residual)
882  {
883  auto & residual = getVector(residualVectorTag());
884 
885  *_residual_ghosted = residual;
887  }
888  // Need to close and update the aux system in case residuals were saved to it.
891  if (hasSaveIn())
893  }
894  catch (MooseException & e)
895  {
896  // The buck stops here, we have already handled the exception by
897  // calling stopSolve(), it is now up to PETSc to return a
898  // "diverged" reason during the next solve.
899  }
900 
901  // not supposed to do anything on matrix
903 
905 }
906 
907 void
908 NonlinearSystemBase::computeResidualAndJacobianTags(const std::set<TagID> & vector_tags,
909  const std::set<TagID> & matrix_tags)
910 {
911  const bool required_residual =
912  vector_tags.find(residualVectorTag()) == vector_tags.end() ? false : true;
913 
914  try
915  {
916  zeroTaggedVectors(vector_tags);
917  computeResidualAndJacobianInternal(vector_tags, matrix_tags);
918  closeTaggedVectors(vector_tags);
919  closeTaggedMatrices(matrix_tags);
920 
921  if (required_residual)
922  {
923  auto & residual = getVector(residualVectorTag());
924  if (!_time_integrators.empty())
925  {
926  for (auto & ti : _time_integrators)
927  ti->postResidual(residual);
928  }
929  else
930  residual += *_Re_non_time;
931  residual.close();
932  }
933 
935  closeTaggedVectors(vector_tags);
936  closeTaggedMatrices(matrix_tags);
937  }
938  catch (MooseException & e)
939  {
940  // The buck stops here, we have already handled the exception by
941  // calling stopSolve(), it is now up to PETSc to return a
942  // "diverged" reason during the next solve.
943  }
944 }
945 
946 void
948 {
949  for (auto & ti : _time_integrators)
950  ti->preSolve();
951  if (_predictor.get())
952  _predictor->timestepSetup();
953 }
954 
955 void
957 {
959 
960  NumericVector<Number> & initial_solution(solution());
961  if (_predictor.get())
962  {
963  if (_predictor->shouldApply())
964  {
965  TIME_SECTION("applyPredictor", 2, "Applying Predictor");
966 
967  _predictor->apply(initial_solution);
968  _fe_problem.predictorCleanup(initial_solution);
969  }
970  else
971  _console << " Skipping predictor this step" << std::endl;
972  }
973 
974  // do nodal BC
975  {
976  TIME_SECTION("initialBCs", 2, "Applying BCs To Initial Condition");
977 
979  for (const auto & bnode : bnd_nodes)
980  {
981  BoundaryID boundary_id = bnode->_bnd_id;
982  Node * node = bnode->_node;
983 
984  if (node->processor_id() == processor_id())
985  {
986  bool has_preset_nodal_bcs = _preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
987  bool has_ad_preset_nodal_bcs = _ad_preset_nodal_bcs.hasActiveBoundaryObjects(boundary_id);
988 
989  // reinit variables in nodes
990  if (has_preset_nodal_bcs || has_ad_preset_nodal_bcs)
991  _fe_problem.reinitNodeFace(node, boundary_id, 0);
992 
993  if (has_preset_nodal_bcs)
994  {
995  const auto & preset_bcs = _preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
996  for (const auto & preset_bc : preset_bcs)
997  preset_bc->computeValue(initial_solution);
998  }
999  if (has_ad_preset_nodal_bcs)
1000  {
1001  const auto & preset_bcs_res = _ad_preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
1002  for (const auto & preset_bc : preset_bcs_res)
1003  preset_bc->computeValue(initial_solution);
1004  }
1005  }
1006  }
1007  }
1008 
1009 #ifdef MOOSE_KOKKOS_ENABLED
1012 #endif
1013 
1014  _sys.solution->close();
1015  update();
1016 
1017  // Set constraint secondary values
1018  setConstraintSecondaryValues(initial_solution, false);
1019 
1021  setConstraintSecondaryValues(initial_solution, true);
1022 }
1023 
1024 void
1025 NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
1026 {
1027  _predictor = predictor;
1028 }
1029 
1030 void
1032 {
1034 
1035  _kernels.subdomainSetup(subdomain, tid);
1036  _nodal_kernels.subdomainSetup(subdomain, tid);
1037  _element_dampers.subdomainSetup(subdomain, tid);
1038  _nodal_dampers.subdomainSetup(subdomain, tid);
1039 }
1040 
1043 {
1044  if (!_Re_time)
1045  {
1047 
1048  // Most applications don't need the expense of ghosting
1050  _Re_time = &addVector(_Re_time_tag, false, ptype);
1051  }
1052  else if (_need_residual_ghosted && _Re_time->type() == PARALLEL)
1053  {
1054  const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
1055 
1056  // If an application changes its mind, the libMesh API lets us
1057  // change the vector.
1058  _Re_time = &system().add_vector(vector_name, false, GHOSTED);
1059  }
1060 
1061  return *_Re_time;
1062 }
1063 
1066 {
1067  if (!_Re_non_time)
1068  {
1070 
1071  // Most applications don't need the expense of ghosting
1073  _Re_non_time = &addVector(_Re_non_time_tag, false, ptype);
1074  }
1076  {
1077  const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
1078 
1079  // If an application changes its mind, the libMesh API lets us
1080  // change the vector.
1081  _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
1082  }
1083 
1084  return *_Re_non_time;
1085 }
1086 
1089 {
1090  mooseDeprecated("Please use getVector()");
1091  switch (tag)
1092  {
1093  case 0:
1094  return getResidualNonTimeVector();
1095 
1096  case 1:
1097  return getResidualTimeVector();
1098 
1099  default:
1100  mooseError("The required residual vector is not available");
1101  }
1102 }
1103 
1104 void
1106 {
1107  THREAD_ID tid = 0; // constraints are going to be done single-threaded
1108  residual.close();
1110  {
1111  const auto & ncs = _constraints.getActiveNodalConstraints();
1112  for (const auto & nc : ncs)
1113  {
1114  std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1115  std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1116 
1117  if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1118  {
1119  _fe_problem.reinitNodes(primary_node_ids, tid);
1120  _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1121  nc->computeResidual(residual);
1122  }
1123  }
1124  _fe_problem.addCachedResidualDirectly(residual, tid);
1125  residual.close();
1126  }
1127 }
1128 
1129 void
1131 {
1132  if (!hasMatrix(systemMatrixTag()))
1133  mooseError(" A system matrix is required");
1134 
1135  auto & jacobian = getMatrix(systemMatrixTag());
1136  THREAD_ID tid = 0; // constraints are going to be done single-threaded
1137 
1139  {
1140  const auto & ncs = _constraints.getActiveNodalConstraints();
1141  for (const auto & nc : ncs)
1142  {
1143  std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
1144  std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
1145 
1146  if ((secondary_node_ids.size() > 0) && (primary_node_ids.size() > 0))
1147  {
1148  _fe_problem.reinitNodes(primary_node_ids, tid);
1149  _fe_problem.reinitNodesNeighbor(secondary_node_ids, tid);
1150  nc->computeJacobian(jacobian);
1151  }
1152  }
1154  }
1155 }
1156 
1157 void
1159  const BoundaryID secondary_boundary,
1160  const PenetrationInfo & info,
1161  const bool displaced)
1162 {
1163  auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1164  : static_cast<SubProblem &>(_fe_problem);
1165 
1166  const Elem * primary_elem = info._elem;
1167  unsigned int primary_side = info._side_num;
1168  std::vector<Point> points;
1169  points.push_back(info._closest_point);
1170 
1171  // *These next steps MUST be done in this order!*
1172  // ADL: This is a Chesterton's fence situation. I don't know which calls exactly the above comment
1173  // is referring to. If I had to guess I would guess just the reinitNodeFace and prepareAssembly
1174  // calls since the former will size the variable's dof indices and then the latter will resize the
1175  // residual/Jacobian based off the variable's cached dof indices size
1176 
1177  // This reinits the variables that exist on the secondary node
1178  _fe_problem.reinitNodeFace(&secondary_node, secondary_boundary, 0);
1179 
1180  // This will set aside residual and jacobian space for the variables that have dofs on
1181  // the secondary node
1183 
1184  _fe_problem.setNeighborSubdomainID(primary_elem, 0);
1185 
1186  //
1187  // Reinit material on undisplaced mesh
1188  //
1189 
1190  const Elem * const undisplaced_primary_elem =
1191  displaced ? _mesh.elemPtr(primary_elem->id()) : primary_elem;
1192  const Point undisplaced_primary_physical_point =
1193  [&points, displaced, primary_elem, undisplaced_primary_elem]()
1194  {
1195  if (displaced)
1196  {
1197  const Point reference_point =
1198  FEMap::inverse_map(primary_elem->dim(), primary_elem, points[0]);
1199  return FEMap::map(primary_elem->dim(), undisplaced_primary_elem, reference_point);
1200  }
1201  else
1202  // If our penetration locator is on the reference mesh, then our undisplaced
1203  // physical point is simply the point coming from the penetration locator
1204  return points[0];
1205  }();
1206 
1208  undisplaced_primary_elem, primary_side, {undisplaced_primary_physical_point}, 0);
1209  // Stateful material properties are only initialized for neighbor material data for internal faces
1210  // for discontinuous Galerkin methods or for conforming interfaces for interface kernels. We don't
1211  // have either of those use cases here where we likely have disconnected meshes
1212  _fe_problem.reinitMaterialsNeighbor(primary_elem->subdomain_id(), 0, /*swap_stateful=*/false);
1213 
1214  // Reinit points for constraint enforcement
1215  if (displaced)
1216  subproblem.reinitNeighborPhys(primary_elem, primary_side, points, 0);
1217 }
1218 
1219 void
1221 {
1222 
1223  if (displaced)
1224  mooseAssert(_fe_problem.getDisplacedProblem(),
1225  "If we're calling this method with displaced = true, then we better well have a "
1226  "displaced problem");
1227  auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1228  : static_cast<SubProblem &>(_fe_problem);
1229  const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1230 
1231  bool constraints_applied = false;
1232 
1233  for (const auto & it : penetration_locators)
1234  {
1235  PenetrationLocator & pen_loc = *(it.second);
1236 
1237  std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1238 
1239  BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1240  BoundaryID primary_boundary = pen_loc._primary_boundary;
1241 
1242  if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1243  {
1244  const auto & constraints =
1245  _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1246  std::unordered_set<unsigned int> needed_mat_props;
1247  for (const auto & constraint : constraints)
1248  {
1249  const auto & mp_deps = constraint->getMatPropDependencies();
1250  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
1251  }
1252  _fe_problem.setActiveMaterialProperties(needed_mat_props, /*tid=*/0);
1253 
1254  for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1255  {
1256  dof_id_type secondary_node_num = secondary_nodes[i];
1257  Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1258 
1259  if (secondary_node.processor_id() == processor_id())
1260  {
1261  if (pen_loc._penetration_info[secondary_node_num])
1262  {
1263  PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1264 
1265  reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1266 
1267  for (const auto & nfc : constraints)
1268  {
1269  if (nfc->isExplicitConstraint())
1270  continue;
1271  // Return if this constraint does not correspond to the primary-secondary pair
1272  // prepared by the outer loops.
1273  // This continue statement is required when, e.g. one secondary surface constrains
1274  // more than one primary surface.
1275  if (nfc->secondaryBoundary() != secondary_boundary ||
1276  nfc->primaryBoundary() != primary_boundary)
1277  continue;
1278 
1279  if (nfc->shouldApply())
1280  {
1281  constraints_applied = true;
1282  nfc->computeSecondaryValue(solution);
1283  }
1284 
1285  if (nfc->hasWritableCoupledVariables())
1286  {
1287  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1288  for (auto * var : nfc->getWritableCoupledVariables())
1289  {
1290  if (var->isNodalDefined())
1291  var->insert(_fe_problem.getAuxiliarySystem().solution());
1292  }
1293  }
1294  }
1295  }
1296  }
1297  }
1298  }
1299  }
1300 
1301  // go over NodeELemConstraints
1302  std::set<dof_id_type> unique_secondary_node_ids;
1303 
1304  for (const auto & secondary_id : _mesh.meshSubdomains())
1305  {
1306  for (const auto & primary_id : _mesh.meshSubdomains())
1307  {
1308  if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1309  {
1310  const auto & constraints =
1311  _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1312 
1313  // get unique set of ids of all nodes on current block
1314  unique_secondary_node_ids.clear();
1315  const MeshBase & meshhelper = _mesh.getMesh();
1316  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1317  meshhelper.active_subdomain_elements_end(secondary_id)))
1318  {
1319  for (auto & n : elem->node_ref_range())
1320  unique_secondary_node_ids.insert(n.id());
1321  }
1322 
1323  for (auto secondary_node_id : unique_secondary_node_ids)
1324  {
1325  Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1326 
1327  // check if secondary node is on current processor
1328  if (secondary_node.processor_id() == processor_id())
1329  {
1330  // This reinits the variables that exist on the secondary node
1331  _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
1332 
1333  // This will set aside residual and jacobian space for the variables that have dofs
1334  // on the secondary node
1336 
1337  for (const auto & nec : constraints)
1338  {
1339  if (nec->shouldApply())
1340  {
1341  constraints_applied = true;
1342  nec->computeSecondaryValue(solution);
1343  }
1344  }
1345  }
1346  }
1347  }
1348  }
1349  }
1350 
1351  // See if constraints were applied anywhere
1352  _communicator.max(constraints_applied);
1353 
1354  if (constraints_applied)
1355  {
1356  solution.close();
1357  update();
1358  }
1359 }
1360 
1361 void
1363 {
1364  // Make sure the residual is in a good state
1365  residual.close();
1366 
1367  if (displaced)
1368  mooseAssert(_fe_problem.getDisplacedProblem(),
1369  "If we're calling this method with displaced = true, then we better well have a "
1370  "displaced problem");
1371  auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1372  : static_cast<SubProblem &>(_fe_problem);
1373  const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1374 
1375  bool constraints_applied;
1376  bool residual_has_inserted_values = false;
1378  constraints_applied = false;
1379  for (const auto & it : penetration_locators)
1380  {
1382  {
1383  // Reset the constraint_applied flag before each new constraint, as they need to be
1384  // assembled separately
1385  constraints_applied = false;
1386  }
1387  PenetrationLocator & pen_loc = *(it.second);
1388 
1389  std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1390 
1391  BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1392  BoundaryID primary_boundary = pen_loc._primary_boundary;
1393 
1394  bool has_writable_variables(false);
1395 
1396  if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
1397  {
1398  const auto & constraints =
1399  _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
1400 
1401  for (unsigned int i = 0; i < secondary_nodes.size(); i++)
1402  {
1403  dof_id_type secondary_node_num = secondary_nodes[i];
1404  Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1405 
1406  if (secondary_node.processor_id() == processor_id())
1407  {
1408  if (pen_loc._penetration_info[secondary_node_num])
1409  {
1410  PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
1411 
1412  reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
1413 
1414  for (const auto & nfc : constraints)
1415  {
1416  // Return if this constraint does not correspond to the primary-secondary pair
1417  // prepared by the outer loops.
1418  // This continue statement is required when, e.g. one secondary surface constrains
1419  // more than one primary surface.
1420  if (nfc->secondaryBoundary() != secondary_boundary ||
1421  nfc->primaryBoundary() != primary_boundary)
1422  continue;
1423 
1424  if (nfc->shouldApply())
1425  {
1426  constraints_applied = true;
1427  nfc->computeResidual();
1428 
1429  if (nfc->overwriteSecondaryResidual())
1430  {
1431  // The below will actually overwrite the residual for every single dof that
1432  // lives on the node. We definitely don't want to do that!
1433  // _fe_problem.setResidual(residual, 0);
1434 
1435  const auto & secondary_var = nfc->variable();
1436  const auto & secondary_dofs = secondary_var.dofIndices();
1437  mooseAssert(secondary_dofs.size() == secondary_var.count(),
1438  "We are on a node so there should only be one dof per variable (for "
1439  "an ArrayVariable we should have a number of dofs equal to the "
1440  "number of components");
1441 
1442  // Assume that if the user is overwriting the secondary residual, then they are
1443  // supplying residuals that do not correspond to their other physics
1444  // (e.g. Kernels), hence we should not apply a scalingFactor that is normally
1445  // based on the order of their other physics (e.g. Kernels)
1446  std::vector<Number> values = {nfc->secondaryResidual()};
1447  residual.insert(values, secondary_dofs);
1448  residual_has_inserted_values = true;
1449  }
1450  else
1453  }
1454  if (nfc->hasWritableCoupledVariables())
1455  {
1456  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1457  has_writable_variables = true;
1458  for (auto * var : nfc->getWritableCoupledVariables())
1459  {
1460  if (var->isNodalDefined())
1461  var->insert(_fe_problem.getAuxiliarySystem().solution());
1462  }
1463  }
1464  }
1465  }
1466  }
1467  }
1468  }
1469  _communicator.max(has_writable_variables);
1470 
1471  if (has_writable_variables)
1472  {
1473  // Explicit contact dynamic constraints write to auxiliary variables and update the old
1474  // displacement solution on the constraint boundaries. Close solutions and update system
1475  // accordingly.
1478  solutionOld().close();
1479  }
1480 
1482  {
1483  // Make sure that secondary contribution to primary are assembled, and ghosts have been
1484  // exchanged, as current primaries might become secondaries on next iteration and will need to
1485  // contribute their former secondaries' contributions to the future primaries. See if
1486  // constraints were applied anywhere
1487  _communicator.max(constraints_applied);
1488 
1489  if (constraints_applied)
1490  {
1491  // If any of the above constraints inserted values in the residual, it needs to be
1492  // assembled before adding the cached residuals below.
1493  _communicator.max(residual_has_inserted_values);
1494  if (residual_has_inserted_values)
1495  {
1496  residual.close();
1497  residual_has_inserted_values = false;
1498  }
1500  residual.close();
1501 
1503  *_residual_ghosted = residual;
1504  }
1505  }
1506  }
1508  {
1509  _communicator.max(constraints_applied);
1510 
1511  if (constraints_applied)
1512  {
1513  // If any of the above constraints inserted values in the residual, it needs to be assembled
1514  // before adding the cached residuals below.
1515  _communicator.max(residual_has_inserted_values);
1516  if (residual_has_inserted_values)
1517  residual.close();
1518 
1520  residual.close();
1521 
1523  *_residual_ghosted = residual;
1524  }
1525  }
1526 
1527  // go over element-element constraint interface
1528  THREAD_ID tid = 0;
1529  const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
1530  for (const auto & it : element_pair_locators)
1531  {
1532  ElementPairLocator & elem_pair_loc = *(it.second);
1533 
1534  if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
1535  {
1536  // ElemElemConstraint objects
1537  const auto & element_constraints =
1538  _constraints.getActiveElemElemConstraints(it.first, displaced);
1539 
1540  // go over pair elements
1541  const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
1542  elem_pair_loc.getElemPairs();
1543  for (const auto & pr : elem_pairs)
1544  {
1545  const Elem * elem1 = pr.first;
1546  const Elem * elem2 = pr.second;
1547 
1548  if (elem1->processor_id() != processor_id())
1549  continue;
1550 
1551  const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
1552 
1553  // for each element process constraints on the
1554  for (const auto & ec : element_constraints)
1555  {
1556  _fe_problem.setCurrentSubdomainID(elem1, tid);
1557  subproblem.reinitElemPhys(elem1, info._elem1_constraint_q_point, tid);
1558  _fe_problem.setNeighborSubdomainID(elem2, tid);
1559  subproblem.reinitNeighborPhys(elem2, info._elem2_constraint_q_point, tid);
1560 
1561  ec->prepareShapes(ec->variable().number());
1562  ec->prepareNeighborShapes(ec->variable().number());
1563 
1564  ec->reinit(info);
1565  ec->computeResidual();
1568  }
1570  }
1571  }
1572  }
1573 
1574  // go over NodeElemConstraints
1575  std::set<dof_id_type> unique_secondary_node_ids;
1576 
1577  constraints_applied = false;
1578  residual_has_inserted_values = false;
1579  bool has_writable_variables = false;
1580  for (const auto & secondary_id : _mesh.meshSubdomains())
1581  {
1582  for (const auto & primary_id : _mesh.meshSubdomains())
1583  {
1584  if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
1585  {
1586  const auto & constraints =
1587  _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
1588 
1589  // get unique set of ids of all nodes on current block
1590  unique_secondary_node_ids.clear();
1591  const MeshBase & meshhelper = _mesh.getMesh();
1592  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
1593  meshhelper.active_subdomain_elements_end(secondary_id)))
1594  {
1595  for (auto & n : elem->node_ref_range())
1596  unique_secondary_node_ids.insert(n.id());
1597  }
1598 
1599  for (auto secondary_node_id : unique_secondary_node_ids)
1600  {
1601  Node & secondary_node = _mesh.nodeRef(secondary_node_id);
1602  // check if secondary node is on current processor
1603  if (secondary_node.processor_id() == processor_id())
1604  {
1605  // This reinits the variables that exist on the secondary node
1606  _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
1607 
1608  // This will set aside residual and jacobian space for the variables that have dofs
1609  // on the secondary node
1611 
1612  for (const auto & nec : constraints)
1613  {
1614  if (nec->shouldApply())
1615  {
1616  constraints_applied = true;
1617  nec->computeResidual();
1618 
1619  if (nec->overwriteSecondaryResidual())
1620  {
1621  _fe_problem.setResidual(residual, 0);
1622  residual_has_inserted_values = true;
1623  }
1624  else
1627  }
1628  if (nec->hasWritableCoupledVariables())
1629  {
1630  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1631  has_writable_variables = true;
1632  for (auto * var : nec->getWritableCoupledVariables())
1633  {
1634  if (var->isNodalDefined())
1635  var->insert(_fe_problem.getAuxiliarySystem().solution());
1636  }
1637  }
1638  }
1640  }
1641  }
1642  }
1643  }
1644  }
1645  _communicator.max(constraints_applied);
1646 
1647  if (constraints_applied)
1648  {
1649  // If any of the above constraints inserted values in the residual, it needs to be assembled
1650  // before adding the cached residuals below.
1651  _communicator.max(residual_has_inserted_values);
1652  if (residual_has_inserted_values)
1653  residual.close();
1654 
1656  residual.close();
1657 
1659  *_residual_ghosted = residual;
1660  }
1661  _communicator.max(has_writable_variables);
1662 
1663  if (has_writable_variables)
1664  {
1665  // Explicit contact dynamic constraints write to auxiliary variables and update the old
1666  // displacement solution on the constraint boundaries. Close solutions and update system
1667  // accordingly.
1670  solutionOld().close();
1671  }
1672 
1673  // We may have additional tagged vectors that also need to be accumulated
1675 }
1676 
1677 void
1679 {
1680  // Overwrite results from integrator in case we have explicit dynamics contact constraints
1682  ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
1683  : static_cast<SubProblem &>(_fe_problem);
1684  const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
1685 
1686  for (const auto & it : penetration_locators)
1687  {
1688  PenetrationLocator & pen_loc = *(it.second);
1689 
1690  const auto & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
1691  const BoundaryID secondary_boundary = pen_loc._secondary_boundary;
1692  const BoundaryID primary_boundary = pen_loc._primary_boundary;
1693 
1694  if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, true))
1695  {
1696  const auto & constraints =
1697  _constraints.getActiveNodeFaceConstraints(secondary_boundary, true);
1698  for (const auto i : index_range(secondary_nodes))
1699  {
1700  const auto secondary_node_num = secondary_nodes[i];
1701  const Node & secondary_node = _mesh.nodeRef(secondary_node_num);
1702 
1703  if (secondary_node.processor_id() == processor_id())
1704  if (pen_loc._penetration_info[secondary_node_num])
1705  for (const auto & nfc : constraints)
1706  {
1707  if (!nfc->isExplicitConstraint())
1708  continue;
1709 
1710  // Return if this constraint does not correspond to the primary-secondary pair
1711  // prepared by the outer loops.
1712  // This continue statement is required when, e.g. one secondary surface constrains
1713  // more than one primary surface.
1714  if (nfc->secondaryBoundary() != secondary_boundary ||
1715  nfc->primaryBoundary() != primary_boundary)
1716  continue;
1717 
1718  nfc->overwriteBoundaryVariables(soln, secondary_node);
1719  }
1720  }
1721  }
1722  }
1723  soln.close();
1724 }
1725 
1726 void
1728 {
1729  TIME_SECTION("residualSetup", 3);
1730 
1732 
1733  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1734  {
1735  _kernels.residualSetup(tid);
1738  if (_doing_dg)
1744  }
1751 
1752 #ifdef MOOSE_KOKKOS_ENABLED
1757 #endif
1758 
1759  // Avoid recursion
1760  if (this == &_fe_problem.currentNonlinearSystem())
1763 }
1764 
1765 void
1767 {
1768  parallel_object_only();
1769 
1770  TIME_SECTION("computeResidualInternal", 3);
1771 
1772  residualSetup();
1773 
1774 #ifdef MOOSE_KOKKOS_ENABLED
1776  computeKokkosResidual(tags);
1777 #endif
1778 
1779  const auto vector_tag_data = _fe_problem.getVectorTags(tags);
1780 
1781  // Residual contributions from UOs - for now this is used for ray tracing
1782  // and ray kernels that contribute to the residual (think line sources)
1783  std::vector<UserObject *> uos;
1785  .query()
1786  .condition<AttribSystem>("UserObject")
1787  .condition<AttribExecOns>(EXEC_PRE_KERNELS)
1788  .queryInto(uos);
1789  for (auto & uo : uos)
1790  uo->residualSetup();
1791  for (auto & uo : uos)
1792  {
1793  uo->initialize();
1794  uo->execute();
1795  uo->finalize();
1796  }
1797 
1798  // reinit scalar variables
1799  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1801 
1802  // residual contributions from the domain
1803  PARALLEL_TRY
1804  {
1805  TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
1806 
1808 
1810  Threads::parallel_reduce(elem_range, cr);
1811 
1812  // We pass face information directly to FV residual objects for their evaluation. Consequently
1813  // we must make sure to do separate threaded loops for 1) undisplaced face information objects
1814  // and undisplaced residual objects and 2) displaced face information objects and displaced
1815  // residual objects
1817  if (_fe_problem.haveFV())
1818  {
1820  _fe_problem, this->number(), tags, /*on_displaced=*/false);
1822  Threads::parallel_reduce(faces, fvr);
1823  }
1825  displaced_problem && displaced_problem->haveFV())
1826  {
1828  _fe_problem, this->number(), tags, /*on_displaced=*/true);
1829  FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
1830  displaced_problem->mesh().ownedFaceInfoEnd());
1831  Threads::parallel_reduce(faces, fvr);
1832  }
1833 
1834  unsigned int n_threads = libMesh::n_threads();
1835  for (unsigned int i = 0; i < n_threads;
1836  i++) // Add any cached residuals that might be hanging around
1838  }
1839  PARALLEL_CATCH;
1840 
1841  // residual contributions from the scalar kernels
1842  PARALLEL_TRY
1843  {
1844  // do scalar kernels (not sure how to thread this)
1846  {
1847  TIME_SECTION("ScalarKernels", 3 /*, "Computing ScalarKernels"*/);
1848 
1849  MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
1850  // This code should be refactored once we can do tags for scalar
1851  // kernels
1852  // Should redo this based on Warehouse
1853  if (!tags.size() || tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL))
1854  scalar_kernel_warehouse = &_scalar_kernels;
1855  else if (tags.size() == 1)
1856  scalar_kernel_warehouse =
1857  &(_scalar_kernels.getVectorTagObjectWarehouse(*(tags.begin()), 0));
1858  else
1859  // scalar_kernels is not threading
1860  scalar_kernel_warehouse = &(_scalar_kernels.getVectorTagsObjectWarehouse(tags, 0));
1861 
1862  bool have_scalar_contributions = false;
1863  const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
1864  for (const auto & scalar_kernel : scalars)
1865  {
1866  scalar_kernel->reinit();
1867  const std::vector<dof_id_type> & dof_indices = scalar_kernel->variable().dofIndices();
1868  const DofMap & dof_map = scalar_kernel->variable().dofMap();
1869  const dof_id_type first_dof = dof_map.first_dof();
1870  const dof_id_type end_dof = dof_map.end_dof();
1871  for (dof_id_type dof : dof_indices)
1872  {
1873  if (dof >= first_dof && dof < end_dof)
1874  {
1875  scalar_kernel->computeResidual();
1876  have_scalar_contributions = true;
1877  break;
1878  }
1879  }
1880  }
1881  if (have_scalar_contributions)
1883  }
1884  }
1885  PARALLEL_CATCH;
1886 
1887  // residual contributions from Block NodalKernels
1888  PARALLEL_TRY
1889  {
1891  {
1892  TIME_SECTION("NodalKernels", 3 /*, "Computing NodalKernels"*/);
1893 
1895 
1897 
1898  if (range.begin() != range.end())
1899  {
1900  _fe_problem.reinitNode(*range.begin(), 0);
1901 
1902  Threads::parallel_reduce(range, cnk);
1903 
1904  unsigned int n_threads = libMesh::n_threads();
1905  for (unsigned int i = 0; i < n_threads;
1906  i++) // Add any cached residuals that might be hanging around
1908  }
1909  }
1910  }
1911  PARALLEL_CATCH;
1912 
1914  // We computed the volumetric objects. We can return now before we get into
1915  // any strongly enforced constraint conditions or penalty-type objects
1916  // (DGKernels, IntegratedBCs, InterfaceKernels, Constraints)
1917  return;
1918 
1919  // residual contributions from boundary NodalKernels
1920  PARALLEL_TRY
1921  {
1923  {
1924  TIME_SECTION("NodalKernelBCs", 3 /*, "Computing NodalKernelBCs"*/);
1925 
1927 
1929 
1930  Threads::parallel_reduce(bnd_node_range, cnk);
1931 
1932  unsigned int n_threads = libMesh::n_threads();
1933  for (unsigned int i = 0; i < n_threads;
1934  i++) // Add any cached residuals that might be hanging around
1936  }
1937  }
1938  PARALLEL_CATCH;
1939 
1941 
1942  if (_residual_copy.get())
1943  {
1944  _Re_non_time->close();
1946  }
1947 
1949  {
1950  _Re_non_time->close();
1953  }
1954 
1955  PARALLEL_TRY { computeDiracContributions(tags, false); }
1956  PARALLEL_CATCH;
1957 
1959  {
1960  PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1961  PARALLEL_CATCH;
1962  _Re_non_time->close();
1963  }
1964 
1965  // Add in Residual contributions from other Constraints
1967  {
1968  PARALLEL_TRY
1969  {
1970  // Undisplaced Constraints
1972 
1973  // Displaced Constraints
1976 
1979  }
1980  PARALLEL_CATCH;
1981  _Re_non_time->close();
1982  }
1983 
1984  // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
1985  // counters
1988 }
1989 
1990 void
1992  const std::set<TagID> & matrix_tags)
1993 {
1994  TIME_SECTION("computeResidualAndJacobianInternal", 3);
1995 
1996  // Make matrix ready to use
1998 
1999  for (auto tag : matrix_tags)
2000  {
2001  if (!hasMatrix(tag))
2002  continue;
2003 
2004  auto & jacobian = getMatrix(tag);
2005  // Necessary for speed
2006  if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2007  {
2008  LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2009  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2010  PETSC_TRUE));
2012  LibmeshPetscCall(
2013  MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2015  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2016  MAT_IGNORE_ZERO_ENTRIES,
2017  PETSC_TRUE));
2018  }
2019  }
2020 
2021  residualSetup();
2022 
2023  // Residual contributions from UOs - for now this is used for ray tracing
2024  // and ray kernels that contribute to the residual (think line sources)
2025  std::vector<UserObject *> uos;
2027  .query()
2028  .condition<AttribSystem>("UserObject")
2029  .condition<AttribExecOns>(EXEC_PRE_KERNELS)
2030  .queryInto(uos);
2031  for (auto & uo : uos)
2032  uo->residualSetup();
2033  for (auto & uo : uos)
2034  {
2035  uo->initialize();
2036  uo->execute();
2037  uo->finalize();
2038  }
2039 
2040  // reinit scalar variables
2041  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2043 
2044  // residual contributions from the domain
2045  PARALLEL_TRY
2046  {
2047  TIME_SECTION("Kernels", 3 /*, "Computing Kernels"*/);
2048 
2050 
2051  ComputeResidualAndJacobianThread crj(_fe_problem, vector_tags, matrix_tags);
2052  Threads::parallel_reduce(elem_range, crj);
2053 
2055  if (_fe_problem.haveFV())
2056  {
2058  _fe_problem, this->number(), vector_tags, matrix_tags, /*on_displaced=*/false);
2060  Threads::parallel_reduce(faces, fvrj);
2061  }
2063  displaced_problem && displaced_problem->haveFV())
2064  {
2066  _fe_problem, this->number(), vector_tags, matrix_tags, /*on_displaced=*/true);
2067  FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
2068  displaced_problem->mesh().ownedFaceInfoEnd());
2069  Threads::parallel_reduce(faces, fvr);
2070  }
2071 
2073 
2074  unsigned int n_threads = libMesh::n_threads();
2075  for (unsigned int i = 0; i < n_threads;
2076  i++) // Add any cached residuals that might be hanging around
2077  {
2080  }
2081  }
2082  PARALLEL_CATCH;
2083 }
2084 
2085 void
2087 {
2088  _nl_vector_tags.clear();
2089 
2090  const auto & residual_vector_tags = _fe_problem.getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
2091  for (const auto & residual_vector_tag : residual_vector_tags)
2092  _nl_vector_tags.insert(residual_vector_tag._id);
2093 
2095  computeNodalBCs(residual, _nl_vector_tags);
2097 }
2098 
2099 void
2100 NonlinearSystemBase::computeNodalBCs(NumericVector<Number> & residual, const std::set<TagID> & tags)
2101 {
2103 
2104  computeNodalBCs(tags);
2105 
2107 }
2108 
2109 void
2110 NonlinearSystemBase::computeNodalBCs(const std::set<TagID> & tags)
2111 {
2112  // We need to close the diag_save_in variables on the aux system before NodalBCBases clear the
2113  // dofs on boundary nodes
2114  if (_has_save_in)
2116 
2117  // Select nodal kernels
2118  MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2119 
2120  if (tags.size() == _fe_problem.numVectorTags(Moose::VECTOR_TAG_RESIDUAL) || !tags.size())
2121  nbc_warehouse = &_nodal_bcs;
2122  else if (tags.size() == 1)
2123  nbc_warehouse = &(_nodal_bcs.getVectorTagObjectWarehouse(*(tags.begin()), 0));
2124  else
2125  nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0));
2126 
2127  // Return early if there is no nodal kernel
2128  if (!nbc_warehouse->size())
2129  return;
2130 
2131  PARALLEL_TRY
2132  {
2134 
2135  if (!bnd_nodes.empty())
2136  {
2137  TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2138 
2139  for (const auto & bnode : bnd_nodes)
2140  {
2141  BoundaryID boundary_id = bnode->_bnd_id;
2142  Node * node = bnode->_node;
2143 
2144  if (node->processor_id() == processor_id() &&
2145  nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
2146  {
2147  // reinit variables in nodes
2148  _fe_problem.reinitNodeFace(node, boundary_id, 0);
2149 
2150  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2151  for (const auto & nbc : bcs)
2152  if (nbc->shouldApply())
2153  nbc->computeResidual();
2154  }
2155  }
2156  }
2157  }
2158  PARALLEL_CATCH;
2159 
2160  if (_Re_time)
2161  _Re_time->close();
2162  _Re_non_time->close();
2163 }
2164 
2165 void
2167 {
2168  PARALLEL_TRY
2169  {
2171 
2172  if (!bnd_nodes.empty())
2173  {
2174  TIME_SECTION("NodalBCs", 3 /*, "Computing NodalBCs"*/);
2175 
2176  for (const auto & bnode : bnd_nodes)
2177  {
2178  BoundaryID boundary_id = bnode->_bnd_id;
2179  Node * node = bnode->_node;
2180 
2181  if (node->processor_id() == processor_id())
2182  {
2183  // reinit variables in nodes
2184  _fe_problem.reinitNodeFace(node, boundary_id, 0);
2185  if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
2186  {
2187  const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
2188  for (const auto & nbc : bcs)
2189  if (nbc->shouldApply())
2190  nbc->computeResidualAndJacobian();
2191  }
2192  }
2193  }
2194  }
2195  }
2196  PARALLEL_CATCH;
2197 
2198  // Set the cached NodalBCBase values in the Jacobian matrix
2200 }
2201 
2202 void
2203 NonlinearSystemBase::getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs)
2204 {
2205  const Node & node = _mesh.nodeRef(node_id);
2206  unsigned int s = number();
2207  if (node.has_dofs(s))
2208  {
2209  for (unsigned int v = 0; v < nVariables(); v++)
2210  for (unsigned int c = 0; c < node.n_comp(s, v); c++)
2211  dofs.push_back(node.dof_number(s, v, c));
2212  }
2213 }
2214 
2215 void
2217  GeometricSearchData & geom_search_data,
2218  std::unordered_map<dof_id_type, std::vector<dof_id_type>> & graph)
2219 {
2220  const auto & node_to_elem_map = _mesh.nodeToElemMap();
2221  const auto & nearest_node_locators = geom_search_data._nearest_node_locators;
2222  for (const auto & it : nearest_node_locators)
2223  {
2224  std::vector<dof_id_type> & secondary_nodes = it.second->_secondary_nodes;
2225 
2226  for (const auto & secondary_node : secondary_nodes)
2227  {
2228  std::set<dof_id_type> unique_secondary_indices;
2229  std::set<dof_id_type> unique_primary_indices;
2230 
2231  auto node_to_elem_pair = node_to_elem_map.find(secondary_node);
2232  if (node_to_elem_pair != node_to_elem_map.end())
2233  {
2234  const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
2235 
2236  // Get the dof indices from each elem connected to the node
2237  for (const auto & cur_elem : elems)
2238  {
2239  std::vector<dof_id_type> dof_indices;
2240  dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2241 
2242  for (const auto & dof : dof_indices)
2243  unique_secondary_indices.insert(dof);
2244  }
2245  }
2246 
2247  std::vector<dof_id_type> primary_nodes = it.second->_neighbor_nodes[secondary_node];
2248 
2249  for (const auto & primary_node : primary_nodes)
2250  {
2251  auto primary_node_to_elem_pair = node_to_elem_map.find(primary_node);
2252  mooseAssert(primary_node_to_elem_pair != node_to_elem_map.end(),
2253  "Missing entry in node to elem map");
2254  const std::vector<dof_id_type> & primary_node_elems = primary_node_to_elem_pair->second;
2255 
2256  // Get the dof indices from each elem connected to the node
2257  for (const auto & cur_elem : primary_node_elems)
2258  {
2259  std::vector<dof_id_type> dof_indices;
2260  dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
2261 
2262  for (const auto & dof : dof_indices)
2263  unique_primary_indices.insert(dof);
2264  }
2265  }
2266 
2267  for (const auto & secondary_id : unique_secondary_indices)
2268  for (const auto & primary_id : unique_primary_indices)
2269  {
2270  graph[secondary_id].push_back(primary_id);
2271  graph[primary_id].push_back(secondary_id);
2272  }
2273  }
2274  }
2275 
2276  // handle node-to-node constraints
2277  const auto & ncs = _constraints.getActiveNodalConstraints();
2278  for (const auto & nc : ncs)
2279  {
2280  std::vector<dof_id_type> primary_dofs;
2281  std::vector<dof_id_type> & primary_node_ids = nc->getPrimaryNodeId();
2282  for (const auto & node_id : primary_node_ids)
2283  {
2284  Node * node = _mesh.queryNodePtr(node_id);
2285  if (node && node->processor_id() == this->processor_id())
2286  {
2287  getNodeDofs(node_id, primary_dofs);
2288  }
2289  }
2290 
2291  _communicator.allgather(primary_dofs);
2292 
2293  std::vector<dof_id_type> secondary_dofs;
2294  std::vector<dof_id_type> & secondary_node_ids = nc->getSecondaryNodeId();
2295  for (const auto & node_id : secondary_node_ids)
2296  {
2297  Node * node = _mesh.queryNodePtr(node_id);
2298  if (node && node->processor_id() == this->processor_id())
2299  {
2300  getNodeDofs(node_id, secondary_dofs);
2301  }
2302  }
2303 
2304  _communicator.allgather(secondary_dofs);
2305 
2306  for (const auto & primary_id : primary_dofs)
2307  for (const auto & secondary_id : secondary_dofs)
2308  {
2309  graph[primary_id].push_back(secondary_id);
2310  graph[secondary_id].push_back(primary_id);
2311  }
2312  }
2313 
2314  // Make every entry sorted and unique
2315  for (auto & it : graph)
2316  {
2317  std::vector<dof_id_type> & row = it.second;
2318  std::sort(row.begin(), row.end());
2319  std::vector<dof_id_type>::iterator uit = std::unique(row.begin(), row.end());
2320  row.resize(uit - row.begin());
2321  }
2322 }
2323 
2324 void
2326 {
2327  if (!hasMatrix(systemMatrixTag()))
2328  mooseError("Need a system matrix ");
2329 
2330  // At this point, have no idea how to make
2331  // this work with tag system
2332  auto & jacobian = getMatrix(systemMatrixTag());
2333 
2334  std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
2335 
2336  findImplicitGeometricCouplingEntries(geom_search_data, graph);
2337 
2338  for (const auto & it : graph)
2339  {
2340  dof_id_type dof = it.first;
2341  const auto & row = it.second;
2342 
2343  for (const auto & coupled_dof : row)
2344  jacobian.add(dof, coupled_dof, 0);
2345  }
2346 }
2347 
2348 void
2350  bool displaced)
2351 {
2352  if (!hasMatrix(systemMatrixTag()))
2353  mooseError("A system matrix is required");
2354 
2355  auto & jacobian = getMatrix(systemMatrixTag());
2356 
2358  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2359  MAT_NEW_NONZERO_ALLOCATION_ERR,
2360  PETSC_FALSE));
2362  LibmeshPetscCall(MatSetOption(
2363  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
2364 
2365  std::vector<numeric_index_type> zero_rows;
2366 
2367  if (displaced)
2368  mooseAssert(_fe_problem.getDisplacedProblem(),
2369  "If we're calling this method with displaced = true, then we better well have a "
2370  "displaced problem");
2371  auto & subproblem = displaced ? static_cast<SubProblem &>(*_fe_problem.getDisplacedProblem())
2372  : static_cast<SubProblem &>(_fe_problem);
2373  const auto & penetration_locators = subproblem.geomSearchData()._penetration_locators;
2374 
2375  bool constraints_applied;
2377  constraints_applied = false;
2378  for (const auto & it : penetration_locators)
2379  {
2381  {
2382  // Reset the constraint_applied flag before each new constraint, as they need to be
2383  // assembled separately
2384  constraints_applied = false;
2385  }
2386  PenetrationLocator & pen_loc = *(it.second);
2387 
2388  std::vector<dof_id_type> & secondary_nodes = pen_loc._nearest_node._secondary_nodes;
2389 
2390  BoundaryID secondary_boundary = pen_loc._secondary_boundary;
2391  BoundaryID primary_boundary = pen_loc._primary_boundary;
2392 
2393  zero_rows.clear();
2394  if (_constraints.hasActiveNodeFaceConstraints(secondary_boundary, displaced))
2395  {
2396  const auto & constraints =
2397  _constraints.getActiveNodeFaceConstraints(secondary_boundary, displaced);
2398 
2399  for (const auto & secondary_node_num : secondary_nodes)
2400  {
2401  Node & secondary_node = _mesh.nodeRef(secondary_node_num);
2402 
2403  if (secondary_node.processor_id() == processor_id())
2404  {
2405  if (pen_loc._penetration_info[secondary_node_num])
2406  {
2407  PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
2408 
2409  reinitNodeFace(secondary_node, secondary_boundary, info, displaced);
2411 
2412  for (const auto & nfc : constraints)
2413  {
2414  if (nfc->isExplicitConstraint())
2415  continue;
2416  // Return if this constraint does not correspond to the primary-secondary pair
2417  // prepared by the outer loops.
2418  // This continue statement is required when, e.g. one secondary surface constrains
2419  // more than one primary surface.
2420  if (nfc->secondaryBoundary() != secondary_boundary ||
2421  nfc->primaryBoundary() != primary_boundary)
2422  continue;
2423 
2424  nfc->_jacobian = &jacobian_to_view;
2425 
2426  if (nfc->shouldApply())
2427  {
2428  constraints_applied = true;
2429 
2430  nfc->prepareShapes(nfc->variable().number());
2431  nfc->prepareNeighborShapes(nfc->variable().number());
2432 
2433  nfc->computeJacobian();
2434 
2435  if (nfc->overwriteSecondaryJacobian())
2436  {
2437  // Add this variable's dof's row to be zeroed
2438  zero_rows.push_back(nfc->variable().nodalDofIndex());
2439  }
2440 
2441  std::vector<dof_id_type> secondary_dofs(1, nfc->variable().nodalDofIndex());
2442 
2443  // Assume that if the user is overwriting the secondary Jacobian, then they are
2444  // supplying Jacobians that do not correspond to their other physics
2445  // (e.g. Kernels), hence we should not apply a scalingFactor that is normally
2446  // based on the order of their other physics (e.g. Kernels)
2447  Real scaling_factor =
2448  nfc->overwriteSecondaryJacobian() ? 1. : nfc->variable().scalingFactor();
2449 
2450  // Cache the jacobian block for the secondary side
2451  nfc->addJacobian(_fe_problem.assembly(0, number()),
2452  nfc->_Kee,
2453  secondary_dofs,
2454  nfc->_connected_dof_indices,
2455  scaling_factor);
2456 
2457  // Cache Ken, Kne, Knn
2458  if (nfc->addCouplingEntriesToJacobian())
2459  {
2460  // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2461  // factor when we're overwriting secondary stuff)
2462  nfc->addJacobian(_fe_problem.assembly(0, number()),
2463  nfc->_Ken,
2464  secondary_dofs,
2465  nfc->primaryVariable().dofIndicesNeighbor(),
2466  scaling_factor);
2467 
2468  // Use _connected_dof_indices to get all the correct columns
2469  nfc->addJacobian(_fe_problem.assembly(0, number()),
2470  nfc->_Kne,
2471  nfc->primaryVariable().dofIndicesNeighbor(),
2472  nfc->_connected_dof_indices,
2473  nfc->primaryVariable().scalingFactor());
2474 
2475  // We've handled Ken and Kne, finally handle Knn
2477  }
2478 
2479  // Do the off-diagonals next
2480  const std::vector<MooseVariableFEBase *> coupled_vars = nfc->getCoupledMooseVars();
2481  for (const auto & jvar : coupled_vars)
2482  {
2483  // Only compute jacobians for nonlinear variables
2484  if (jvar->kind() != Moose::VAR_SOLVER)
2485  continue;
2486 
2487  // Only compute Jacobian entries if this coupling is being used by the
2488  // preconditioner
2489  if (nfc->variable().number() == jvar->number() ||
2491  nfc->variable().number(), jvar->number(), this->number()))
2492  continue;
2493 
2494  // Need to zero out the matrices first
2496 
2497  nfc->prepareShapes(nfc->variable().number());
2498  nfc->prepareNeighborShapes(jvar->number());
2499 
2500  nfc->computeOffDiagJacobian(jvar->number());
2501 
2502  // Cache the jacobian block for the secondary side
2503  nfc->addJacobian(_fe_problem.assembly(0, number()),
2504  nfc->_Kee,
2505  secondary_dofs,
2506  nfc->_connected_dof_indices,
2507  scaling_factor);
2508 
2509  // Cache Ken, Kne, Knn
2510  if (nfc->addCouplingEntriesToJacobian())
2511  {
2512  // Make sure we use a proper scaling factor (e.g. don't use an interior scaling
2513  // factor when we're overwriting secondary stuff)
2514  nfc->addJacobian(_fe_problem.assembly(0, number()),
2515  nfc->_Ken,
2516  secondary_dofs,
2517  jvar->dofIndicesNeighbor(),
2518  scaling_factor);
2519 
2520  // Use _connected_dof_indices to get all the correct columns
2521  nfc->addJacobian(_fe_problem.assembly(0, number()),
2522  nfc->_Kne,
2523  nfc->variable().dofIndicesNeighbor(),
2524  nfc->_connected_dof_indices,
2525  nfc->variable().scalingFactor());
2526 
2527  // We've handled Ken and Kne, finally handle Knn
2529  }
2530  }
2531  }
2532  }
2533  }
2534  }
2535  }
2536  }
2538  {
2539  // See if constraints were applied anywhere
2540  _communicator.max(constraints_applied);
2541 
2542  if (constraints_applied)
2543  {
2544  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2545  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2546  PETSC_TRUE));
2547 
2548  jacobian.close();
2549  jacobian.zero_rows(zero_rows, 0.0);
2550  jacobian.close();
2552  jacobian.close();
2553  }
2554  }
2555  }
2557  {
2558  // See if constraints were applied anywhere
2559  _communicator.max(constraints_applied);
2560 
2561  if (constraints_applied)
2562  {
2563  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2564  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2565  PETSC_TRUE));
2566 
2567  jacobian.close();
2568  jacobian.zero_rows(zero_rows, 0.0);
2569  jacobian.close();
2571  jacobian.close();
2572  }
2573  }
2574 
2575  THREAD_ID tid = 0;
2576  // go over element-element constraint interface
2577  const auto & element_pair_locators = subproblem.geomSearchData()._element_pair_locators;
2578  for (const auto & it : element_pair_locators)
2579  {
2580  ElementPairLocator & elem_pair_loc = *(it.second);
2581 
2582  if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
2583  {
2584  // ElemElemConstraint objects
2585  const auto & element_constraints =
2586  _constraints.getActiveElemElemConstraints(it.first, displaced);
2587 
2588  // go over pair elements
2589  const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
2590  elem_pair_loc.getElemPairs();
2591  for (const auto & pr : elem_pairs)
2592  {
2593  const Elem * elem1 = pr.first;
2594  const Elem * elem2 = pr.second;
2595 
2596  if (elem1->processor_id() != processor_id())
2597  continue;
2598 
2599  const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
2600 
2601  // for each element process constraints on the
2602  for (const auto & ec : element_constraints)
2603  {
2604  _fe_problem.setCurrentSubdomainID(elem1, tid);
2605  subproblem.reinitElemPhys(elem1, info._elem1_constraint_q_point, tid);
2606  _fe_problem.setNeighborSubdomainID(elem2, tid);
2607  subproblem.reinitNeighborPhys(elem2, info._elem2_constraint_q_point, tid);
2608 
2609  ec->prepareShapes(ec->variable().number());
2610  ec->prepareNeighborShapes(ec->variable().number());
2611 
2612  ec->reinit(info);
2613  ec->computeJacobian();
2616  }
2618  }
2619  }
2620  }
2621 
2622  // go over NodeElemConstraints
2623  std::set<dof_id_type> unique_secondary_node_ids;
2624  constraints_applied = false;
2625  for (const auto & secondary_id : _mesh.meshSubdomains())
2626  {
2627  for (const auto & primary_id : _mesh.meshSubdomains())
2628  {
2629  if (_constraints.hasActiveNodeElemConstraints(secondary_id, primary_id, displaced))
2630  {
2631  const auto & constraints =
2632  _constraints.getActiveNodeElemConstraints(secondary_id, primary_id, displaced);
2633 
2634  // get unique set of ids of all nodes on current block
2635  unique_secondary_node_ids.clear();
2636  const MeshBase & meshhelper = _mesh.getMesh();
2637  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(secondary_id),
2638  meshhelper.active_subdomain_elements_end(secondary_id)))
2639  {
2640  for (auto & n : elem->node_ref_range())
2641  unique_secondary_node_ids.insert(n.id());
2642  }
2643 
2644  for (auto secondary_node_id : unique_secondary_node_ids)
2645  {
2646  const Node & secondary_node = _mesh.nodeRef(secondary_node_id);
2647  // check if secondary node is on current processor
2648  if (secondary_node.processor_id() == processor_id())
2649  {
2650  // This reinits the variables that exist on the secondary node
2651  _fe_problem.reinitNodeFace(&secondary_node, secondary_id, 0);
2652 
2653  // This will set aside residual and jacobian space for the variables that have dofs
2654  // on the secondary node
2657 
2658  for (const auto & nec : constraints)
2659  {
2660  if (nec->shouldApply())
2661  {
2662  constraints_applied = true;
2663 
2664  nec->_jacobian = &jacobian_to_view;
2665  nec->prepareShapes(nec->variable().number());
2666  nec->prepareNeighborShapes(nec->variable().number());
2667 
2668  nec->computeJacobian();
2669 
2670  if (nec->overwriteSecondaryJacobian())
2671  {
2672  // Add this variable's dof's row to be zeroed
2673  zero_rows.push_back(nec->variable().nodalDofIndex());
2674  }
2675 
2676  std::vector<dof_id_type> secondary_dofs(1, nec->variable().nodalDofIndex());
2677 
2678  // Cache the jacobian block for the secondary side
2679  nec->addJacobian(_fe_problem.assembly(0, number()),
2680  nec->_Kee,
2681  secondary_dofs,
2682  nec->_connected_dof_indices,
2683  nec->variable().scalingFactor());
2684 
2685  // Cache the jacobian block for the primary side
2686  nec->addJacobian(_fe_problem.assembly(0, number()),
2687  nec->_Kne,
2688  nec->primaryVariable().dofIndicesNeighbor(),
2689  nec->_connected_dof_indices,
2690  nec->primaryVariable().scalingFactor());
2691 
2694 
2695  // Do the off-diagonals next
2696  const std::vector<MooseVariableFEBase *> coupled_vars = nec->getCoupledMooseVars();
2697  for (const auto & jvar : coupled_vars)
2698  {
2699  // Only compute jacobians for nonlinear variables
2700  if (jvar->kind() != Moose::VAR_SOLVER)
2701  continue;
2702 
2703  // Only compute Jacobian entries if this coupling is being used by the
2704  // preconditioner
2705  if (nec->variable().number() == jvar->number() ||
2707  nec->variable().number(), jvar->number(), this->number()))
2708  continue;
2709 
2710  // Need to zero out the matrices first
2712 
2713  nec->prepareShapes(nec->variable().number());
2714  nec->prepareNeighborShapes(jvar->number());
2715 
2716  nec->computeOffDiagJacobian(jvar->number());
2717 
2718  // Cache the jacobian block for the secondary side
2719  nec->addJacobian(_fe_problem.assembly(0, number()),
2720  nec->_Kee,
2721  secondary_dofs,
2722  nec->_connected_dof_indices,
2723  nec->variable().scalingFactor());
2724 
2725  // Cache the jacobian block for the primary side
2726  nec->addJacobian(_fe_problem.assembly(0, number()),
2727  nec->_Kne,
2728  nec->variable().dofIndicesNeighbor(),
2729  nec->_connected_dof_indices,
2730  nec->variable().scalingFactor());
2731 
2734  }
2735  }
2736  }
2737  }
2738  }
2739  }
2740  }
2741  }
2742  // See if constraints were applied anywhere
2743  _communicator.max(constraints_applied);
2744 
2745  if (constraints_applied)
2746  {
2747  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2748  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2749  PETSC_TRUE));
2750 
2751  jacobian.close();
2752  jacobian.zero_rows(zero_rows, 0.0);
2753  jacobian.close();
2755  jacobian.close();
2756  }
2757 }
2758 
2759 void
2761 {
2762  MooseObjectWarehouse<ScalarKernelBase> * scalar_kernel_warehouse;
2763 
2764  if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2765  scalar_kernel_warehouse = &_scalar_kernels;
2766  else if (tags.size() == 1)
2767  scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2768  else
2769  scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagsObjectWarehouse(tags, 0));
2770 
2771  // Compute the diagonal block for scalar variables
2772  if (scalar_kernel_warehouse->hasActiveObjects())
2773  {
2774  const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
2775 
2776  _fe_problem.reinitScalars(/*tid=*/0);
2777 
2778  _fe_problem.reinitOffDiagScalars(/*_tid*/ 0);
2779 
2780  bool have_scalar_contributions = false;
2781  for (const auto & kernel : scalars)
2782  {
2783  if (!kernel->computesJacobian())
2784  continue;
2785 
2786  kernel->reinit();
2787  const std::vector<dof_id_type> & dof_indices = kernel->variable().dofIndices();
2788  const DofMap & dof_map = kernel->variable().dofMap();
2789  const dof_id_type first_dof = dof_map.first_dof();
2790  const dof_id_type end_dof = dof_map.end_dof();
2791  for (dof_id_type dof : dof_indices)
2792  {
2793  if (dof >= first_dof && dof < end_dof)
2794  {
2795  kernel->computeJacobian();
2796  _fe_problem.addJacobianOffDiagScalar(kernel->variable().number());
2797  have_scalar_contributions = true;
2798  break;
2799  }
2800  }
2801  }
2802 
2803  if (have_scalar_contributions)
2805  }
2806 }
2807 
2808 void
2810 {
2812 
2813  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2814  {
2815  _kernels.jacobianSetup(tid);
2818  if (_doing_dg)
2824  }
2831 
2832 #ifdef MOOSE_KOKKOS_ENABLED
2837 #endif
2838 
2839  // Avoid recursion
2840  if (this == &_fe_problem.currentNonlinearSystem())
2843 }
2844 
2845 void
2847 {
2848  TIME_SECTION("computeJacobianInternal", 3);
2849 
2851 
2852  // Make matrix ready to use
2854 
2855  for (auto tag : tags)
2856  {
2857  if (!hasMatrix(tag))
2858  continue;
2859 
2860  auto & jacobian = getMatrix(tag);
2861  // Necessary for speed
2862  if (auto petsc_matrix = dynamic_cast<PetscMatrix<Number> *>(&jacobian))
2863  {
2864  LibmeshPetscCall(MatSetOption(petsc_matrix->mat(),
2865  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2866  PETSC_TRUE));
2868  LibmeshPetscCall(
2869  MatSetOption(petsc_matrix->mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
2871  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2872  MAT_IGNORE_ZERO_ENTRIES,
2873  PETSC_TRUE));
2874  }
2875  }
2876 
2877  jacobianSetup();
2878 
2879 #ifdef MOOSE_KOKKOS_ENABLED
2881  computeKokkosJacobian(tags);
2882 #endif
2883 
2884  // Jacobian contributions from UOs - for now this is used for ray tracing
2885  // and ray kernels that contribute to the Jacobian (think line sources)
2886  std::vector<UserObject *> uos;
2888  .query()
2889  .condition<AttribSystem>("UserObject")
2890  .condition<AttribExecOns>(EXEC_PRE_KERNELS)
2891  .queryInto(uos);
2892  for (auto & uo : uos)
2893  uo->jacobianSetup();
2894  for (auto & uo : uos)
2895  {
2896  uo->initialize();
2897  uo->execute();
2898  uo->finalize();
2899  }
2900 
2901  // reinit scalar variables
2902  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2904 
2905  PARALLEL_TRY
2906  {
2907  // We would like to compute ScalarKernels, block NodalKernels, FVFluxKernels, and mortar objects
2908  // up front because we want these included whether we are computing an ordinary Jacobian or a
2909  // Jacobian for determining variable scaling factors
2911 
2912  // Block restricted Nodal Kernels
2914  {
2917  Threads::parallel_reduce(range, cnkjt);
2918 
2919  unsigned int n_threads = libMesh::n_threads();
2920  for (unsigned int i = 0; i < n_threads;
2921  i++) // Add any cached jacobians that might be hanging around
2923  }
2924 
2926  if (_fe_problem.haveFV())
2927  {
2928  // the same loop works for both residual and jacobians because it keys
2929  // off of FEProblem's _currently_computing_jacobian parameter
2931  _fe_problem, this->number(), tags, /*on_displaced=*/false);
2933  Threads::parallel_reduce(faces, fvj);
2934  }
2936  displaced_problem && displaced_problem->haveFV())
2937  {
2939  _fe_problem, this->number(), tags, /*on_displaced=*/true);
2940  FVRange faces(displaced_problem->mesh().ownedFaceInfoBegin(),
2941  displaced_problem->mesh().ownedFaceInfoEnd());
2942  Threads::parallel_reduce(faces, fvr);
2943  }
2944 
2946 
2947  // Get our element range for looping over
2949 
2951  {
2952  // Only compute Jacobians corresponding to the diagonals of volumetric compute objects
2953  // because this typically gives us a good representation of the physics. NodalBCs and
2954  // Constraints can introduce dramatically different scales (often order unity).
2955  // IntegratedBCs and/or InterfaceKernels may use penalty factors. DGKernels may be ok, but
2956  // they are almost always used in conjunction with Kernels
2958  Threads::parallel_reduce(elem_range, cj);
2959  unsigned int n_threads = libMesh::n_threads();
2960  for (unsigned int i = 0; i < n_threads;
2961  i++) // Add any Jacobian contributions still hanging around
2963 
2964  // Check whether any exceptions were thrown and propagate this information for parallel
2965  // consistency before
2966  // 1) we do parallel communication when closing tagged matrices
2967  // 2) early returning before reaching our PARALLEL_CATCH below
2969 
2970  closeTaggedMatrices(tags);
2971 
2972  return;
2973  }
2974 
2975  switch (_fe_problem.coupling())
2976  {
2977  case Moose::COUPLING_DIAG:
2978  {
2980  Threads::parallel_reduce(elem_range, cj);
2981 
2982  unsigned int n_threads = libMesh::n_threads();
2983  for (unsigned int i = 0; i < n_threads;
2984  i++) // Add any Jacobian contributions still hanging around
2986 
2987  // Boundary restricted Nodal Kernels
2989  {
2992 
2993  Threads::parallel_reduce(bnd_range, cnkjt);
2994  unsigned int n_threads = libMesh::n_threads();
2995  for (unsigned int i = 0; i < n_threads;
2996  i++) // Add any cached jacobians that might be hanging around
2998  }
2999  }
3000  break;
3001 
3002  default:
3004  {
3006  Threads::parallel_reduce(elem_range, cj);
3007  unsigned int n_threads = libMesh::n_threads();
3008 
3009  for (unsigned int i = 0; i < n_threads; i++)
3011 
3012  // Boundary restricted Nodal Kernels
3014  {
3017 
3018  Threads::parallel_reduce(bnd_range, cnkjt);
3019  unsigned int n_threads = libMesh::n_threads();
3020  for (unsigned int i = 0; i < n_threads;
3021  i++) // Add any cached jacobians that might be hanging around
3023  }
3024  }
3025  break;
3026  }
3027 
3028  computeDiracContributions(tags, true);
3029 
3030  static bool first = true;
3031 
3032  // This adds zeroes into geometric coupling entries to ensure they stay in the matrix
3033  if ((_fe_problem.restoreOriginalNonzeroPattern() || first) &&
3035  {
3036  first = false;
3038 
3041  }
3042  }
3043  PARALLEL_CATCH;
3044 
3045  // Have no idea how to have constraints work
3046  // with the tag system
3047  PARALLEL_TRY
3048  {
3049  // Add in Jacobian contributions from other Constraints
3050  if (_fe_problem._has_constraints && tags.count(systemMatrixTag()))
3051  {
3052  // Some constraints need to be able to read values from the Jacobian, which requires that it
3053  // be closed/assembled
3054  auto & system_matrix = getMatrix(systemMatrixTag());
3055 #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
3056  SparseMatrix<Number> * view_jac_ptr;
3057  std::unique_ptr<SparseMatrix<Number>> hash_copy;
3058  if (system_matrix.use_hash_table())
3059  {
3060  hash_copy = libMesh::cast_ref<PetscMatrix<Number> &>(system_matrix).copy_from_hash();
3061  view_jac_ptr = hash_copy.get();
3062  }
3063  else
3064  view_jac_ptr = &system_matrix;
3065  auto & jacobian_to_view = *view_jac_ptr;
3066 #else
3067  auto & jacobian_to_view = system_matrix;
3068 #endif
3069  if (&jacobian_to_view == &system_matrix)
3070  system_matrix.close();
3071 
3072  // Nodal Constraints
3074 
3075  // Undisplaced Constraints
3076  constraintJacobians(jacobian_to_view, false);
3077 
3078  // Displaced Constraints
3080  constraintJacobians(jacobian_to_view, true);
3081  }
3082  }
3083  PARALLEL_CATCH;
3084 
3085  // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
3086  // on boundary nodes
3087  if (_has_diag_save_in)
3089 
3090  PARALLEL_TRY
3091  {
3092  MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
3093  // Select nodal kernels
3094  if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
3095  nbc_warehouse = &_nodal_bcs;
3096  else if (tags.size() == 1)
3097  nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
3098  else
3099  nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0));
3100 
3101  if (nbc_warehouse->hasActiveObjects())
3102  {
3103  // We may be switching from add to set. Moreover, we rely on a call to MatZeroRows to enforce
3104  // the nodal boundary condition constraints which requires that the matrix be truly assembled
3105  // as opposed to just flushed. Consequently we can't do the following despite any desire to
3106  // keep our initial sparsity pattern honored (see https://gitlab.com/petsc/petsc/-/issues/852)
3107  //
3108  // flushTaggedMatrices(tags);
3109  closeTaggedMatrices(tags);
3110 
3111  // Cache the information about which BCs are coupled to which
3112  // variables, so we don't have to figure it out for each node.
3113  std::map<std::string, std::set<unsigned int>> bc_involved_vars;
3114  const std::set<BoundaryID> & all_boundary_ids = _mesh.getBoundaryIDs();
3115  for (const auto & bid : all_boundary_ids)
3116  {
3117  // Get reference to all the NodalBCs for this ID. This is only
3118  // safe if there are NodalBCBases there to be gotten...
3119  if (nbc_warehouse->hasActiveBoundaryObjects(bid))
3120  {
3121  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
3122  for (const auto & bc : bcs)
3123  {
3124  const std::vector<MooseVariableFEBase *> & coupled_moose_vars =
3125  bc->getCoupledMooseVars();
3126 
3127  // Create the set of "involved" MOOSE nonlinear vars, which includes all coupled vars
3128  // and the BC's own variable
3129  std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3130  for (const auto & coupled_var : coupled_moose_vars)
3131  if (coupled_var->kind() == Moose::VAR_SOLVER)
3132  var_set.insert(coupled_var->number());
3133 
3134  var_set.insert(bc->variable().number());
3135  }
3136  }
3137  }
3138 
3139  // reinit scalar variables again. This reinit does not re-fill any of the scalar variable
3140  // solution arrays because that was done above. It only will reorder the derivative
3141  // information for AD calculations to be suitable for NodalBC calculations
3142  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3143  _fe_problem.reinitScalars(tid, true);
3144 
3145  // Get variable coupling list. We do all the NodalBCBase stuff on
3146  // thread 0... The couplingEntries() data structure determines
3147  // which variables are "coupled" as far as the preconditioner is
3148  // concerned, not what variables a boundary condition specifically
3149  // depends on.
3150  auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number());
3151 
3152  // Compute Jacobians for NodalBCBases
3154  for (const auto & bnode : bnd_nodes)
3155  {
3156  BoundaryID boundary_id = bnode->_bnd_id;
3157  Node * node = bnode->_node;
3158 
3159  if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
3160  node->processor_id() == processor_id())
3161  {
3162  _fe_problem.reinitNodeFace(node, boundary_id, 0);
3163 
3164  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
3165  for (const auto & bc : bcs)
3166  {
3167  // Get the set of involved MOOSE vars for this BC
3168  std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
3169 
3170  // Loop over all the variables whose Jacobian blocks are
3171  // actually being computed, call computeOffDiagJacobian()
3172  // for each one which is actually coupled (otherwise the
3173  // value is zero.)
3174  for (const auto & it : coupling_entries)
3175  {
3176  unsigned int ivar = it.first->number(), jvar = it.second->number();
3177 
3178  // We are only going to call computeOffDiagJacobian() if:
3179  // 1.) the BC's variable is ivar
3180  // 2.) jvar is "involved" with the BC (including jvar==ivar), and
3181  // 3.) the BC should apply.
3182  if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
3183  bc->computeOffDiagJacobian(jvar);
3184  }
3185 
3186  const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars();
3187  for (const auto & jvariable : coupled_scalar_vars)
3188  if (hasScalarVariable(jvariable->name()))
3189  bc->computeOffDiagJacobianScalar(jvariable->number());
3190  }
3191  }
3192  } // end loop over boundary nodes
3193 
3194  // Set the cached NodalBCBase values in the Jacobian matrix
3196  }
3197  }
3198  PARALLEL_CATCH;
3199 
3200  closeTaggedMatrices(tags);
3201 
3202  // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
3203  // on boundary nodes
3206 
3207  if (hasDiagSaveIn())
3209 
3210  // Accumulate the occurrence of solution invalid warnings for the current iteration cumulative
3211  // counters
3214 }
3215 
3216 void
3218 {
3219  _nl_matrix_tags.clear();
3220 
3221  auto & tags = _fe_problem.getMatrixTags();
3222 
3223  for (auto & tag : tags)
3224  _nl_matrix_tags.insert(tag.second);
3225 
3226  computeJacobian(jacobian, _nl_matrix_tags);
3227 }
3228 
3229 void
3230 NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian, const std::set<TagID> & tags)
3231 {
3233 
3234  computeJacobianTags(tags);
3235 
3237 }
3238 
3239 void
3240 NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
3241 {
3242  TIME_SECTION("computeJacobianTags", 5);
3243 
3244  FloatingPointExceptionGuard fpe_guard(_app);
3245 
3246  try
3247  {
3249  }
3250  catch (MooseException & e)
3251  {
3252  // The buck stops here, we have already handled the exception by
3253  // calling stopSolve(), it is now up to PETSc to return a
3254  // "diverged" reason during the next solve.
3255  }
3256 }
3257 
3258 void
3260 {
3261  _nl_matrix_tags.clear();
3262 
3263  auto & tags = _fe_problem.getMatrixTags();
3264  for (auto & tag : tags)
3265  _nl_matrix_tags.insert(tag.second);
3266 
3268 }
3269 
3270 void
3272  const std::set<TagID> & tags)
3273 {
3274  TIME_SECTION("computeJacobianBlocks", 3);
3275  FloatingPointExceptionGuard fpe_guard(_app);
3276 
3277  for (unsigned int i = 0; i < blocks.size(); i++)
3278  {
3279  SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3280 
3281  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3282  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
3283  PETSC_TRUE));
3285  LibmeshPetscCall(MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
3286  MAT_NEW_NONZERO_ALLOCATION_ERR,
3287  PETSC_TRUE));
3288 
3289  jacobian.zero();
3290  }
3291 
3292  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
3294 
3295  PARALLEL_TRY
3296  {
3299  Threads::parallel_reduce(elem_range, cjb);
3300  }
3301  PARALLEL_CATCH;
3302 
3303  for (unsigned int i = 0; i < blocks.size(); i++)
3304  blocks[i]->_jacobian.close();
3305 
3306  for (unsigned int i = 0; i < blocks.size(); i++)
3307  {
3308  libMesh::System & precond_system = blocks[i]->_precond_system;
3309  SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
3310 
3311  unsigned int ivar = blocks[i]->_ivar;
3312  unsigned int jvar = blocks[i]->_jvar;
3313 
3314  // Dirichlet BCs
3315  std::vector<numeric_index_type> zero_rows;
3316  PARALLEL_TRY
3317  {
3319  for (const auto & bnode : bnd_nodes)
3320  {
3321  BoundaryID boundary_id = bnode->_bnd_id;
3322  Node * node = bnode->_node;
3323 
3324  if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
3325  {
3326  const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
3327 
3328  if (node->processor_id() == processor_id())
3329  {
3330  _fe_problem.reinitNodeFace(node, boundary_id, 0);
3331 
3332  for (const auto & bc : bcs)
3333  if (bc->variable().number() == ivar && bc->shouldApply())
3334  {
3335  // The first zero is for the variable number... there is only one variable in
3336  // each mini-system The second zero only works with Lagrange elements!
3337  zero_rows.push_back(node->dof_number(precond_system.number(), 0, 0));
3338  }
3339  }
3340  }
3341  }
3342  }
3343  PARALLEL_CATCH;
3344 
3345  jacobian.close();
3346 
3347  // This zeroes the rows corresponding to Dirichlet BCs and puts a 1.0 on the diagonal
3348  if (ivar == jvar)
3349  jacobian.zero_rows(zero_rows, 1.0);
3350  else
3351  jacobian.zero_rows(zero_rows, 0.0);
3352 
3353  jacobian.close();
3354  }
3355 }
3356 
3357 void
3359 {
3366  _kernels.updateActive(tid);
3368 
3369  if (tid == 0)
3370  {
3378 
3379 #ifdef MOOSE_KOKKOS_ENABLED
3385 #endif
3386  }
3387 }
3388 
3389 Real
3391  const NumericVector<Number> & update)
3392 {
3393  // Default to no damping
3394  Real damping = 1.0;
3395  bool has_active_dampers = false;
3396 
3397  try
3398  {
3400  {
3401  PARALLEL_TRY
3402  {
3403  TIME_SECTION("computeDampers", 3, "Computing Dampers");
3404  has_active_dampers = true;
3407  Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicElementRange(), cid);
3408  damping = std::min(cid.damping(), damping);
3409  }
3410  PARALLEL_CATCH;
3411  }
3412 
3414  {
3415  PARALLEL_TRY
3416  {
3417  TIME_SECTION("computeDamping::element", 3, "Computing Element Damping");
3418 
3419  has_active_dampers = true;
3422  Threads::parallel_reduce(_fe_problem.getCurrentAlgebraicNodeRange(), cndt);
3423  damping = std::min(cndt.damping(), damping);
3424  }
3425  PARALLEL_CATCH;
3426  }
3427 
3429  {
3430  PARALLEL_TRY
3431  {
3432  TIME_SECTION("computeDamping::general", 3, "Computing General Damping");
3433 
3434  has_active_dampers = true;
3435  const auto & gdampers = _general_dampers.getActiveObjects();
3436  for (const auto & damper : gdampers)
3437  {
3438  Real gd_damping = damper->computeDamping(solution, update);
3439  try
3440  {
3441  damper->checkMinDamping(gd_damping);
3442  }
3443  catch (MooseException & e)
3444  {
3446  }
3447  damping = std::min(gd_damping, damping);
3448  }
3449  }
3450  PARALLEL_CATCH;
3451  }
3452  }
3453  catch (MooseException & e)
3454  {
3455  // The buck stops here, we have already handled the exception by
3456  // calling stopSolve(), it is now up to PETSc to return a
3457  // "diverged" reason during the next solve.
3458  }
3459 
3460  _communicator.min(damping);
3461 
3462  if (has_active_dampers && damping < 1.0)
3463  _console << " Damping factor: " << damping << std::endl;
3464 
3465  return damping;
3466 }
3467 
3468 void
3469 NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
3470 {
3472 
3473  std::set<const Elem *> dirac_elements;
3474 
3476  {
3477  TIME_SECTION("computeDirac", 3, "Computing DiracKernels");
3478 
3479  // TODO: Need a threading fix... but it's complicated!
3480  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3481  {
3482  const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
3483  for (const auto & dkernel : dkernels)
3484  {
3485  dkernel->clearPoints();
3486  dkernel->addPoints();
3487  }
3488  }
3489 
3490  ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
3491 
3492  _fe_problem.getDiracElements(dirac_elements);
3493 
3494  DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
3495  // TODO: Make Dirac work thread!
3496  // Threads::parallel_reduce(range, cd);
3497 
3498  cd(range);
3499  }
3500 }
3501 
3504 {
3505  if (!_residual_copy.get())
3507 
3508  return *_residual_copy;
3509 }
3510 
3513 {
3514  _need_residual_ghosted = true;
3515  if (!_residual_ghosted)
3516  {
3517  // The first time we realize we need a ghosted residual vector,
3518  // we add it.
3519  _residual_ghosted = &addVector("residual_ghosted", false, GHOSTED);
3520 
3521  // If we've already realized we need time and/or non-time
3522  // residual vectors, but we haven't yet realized they need to be
3523  // ghosted, fix that now.
3524  //
3525  // If an application changes its mind, the libMesh API lets us
3526  // change the vector.
3527  if (_Re_time)
3528  {
3529  const auto vector_name = _subproblem.vectorTagName(_Re_time_tag);
3530  _Re_time = &system().add_vector(vector_name, false, GHOSTED);
3531  }
3532  if (_Re_non_time)
3533  {
3534  const auto vector_name = _subproblem.vectorTagName(_Re_non_time_tag);
3535  _Re_non_time = &system().add_vector(vector_name, false, GHOSTED);
3536  }
3537  }
3538  return *_residual_ghosted;
3539 }
3540 
3541 void
3543  std::vector<dof_id_type> & n_nz,
3544  std::vector<dof_id_type> & n_oz)
3545 {
3547  {
3549 
3550  std::unordered_map<dof_id_type, std::vector<dof_id_type>> graph;
3551 
3553 
3556  graph);
3557 
3558  const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
3559  const dof_id_type end_dof_on_proc = dofMap().end_dof(processor_id());
3560 
3561  // The total number of dofs on and off processor
3562  const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
3563  const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
3564 
3565  for (const auto & git : graph)
3566  {
3567  dof_id_type dof = git.first;
3568  dof_id_type local_dof = dof - first_dof_on_proc;
3569 
3570  if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
3571  continue;
3572 
3573  const auto & row = git.second;
3574 
3575  SparsityPattern::Row & sparsity_row = sparsity[local_dof];
3576 
3577  unsigned int original_row_length = sparsity_row.size();
3578 
3579  sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
3580 
3581  SparsityPattern::sort_row(
3582  sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
3583 
3584  // Fix up nonzero arrays
3585  for (const auto & coupled_dof : row)
3586  {
3587  if (coupled_dof < first_dof_on_proc || coupled_dof >= end_dof_on_proc)
3588  {
3589  if (n_oz[local_dof] < n_dofs_not_on_proc)
3590  n_oz[local_dof]++;
3591  }
3592  else
3593  {
3594  if (n_nz[local_dof] < n_dofs_on_proc)
3595  n_nz[local_dof]++;
3596  }
3597  }
3598  }
3599  }
3600 }
3601 
3602 void
3604 {
3605  *_u_dot = u_dot;
3606 }
3607 
3608 void
3610 {
3611  *_u_dotdot = u_dotdot;
3612 }
3613 
3614 void
3616 {
3617  *_u_dot_old = u_dot_old;
3618 }
3619 
3620 void
3622 {
3623  *_u_dotdot_old = u_dotdot_old;
3624 }
3625 
3626 void
3627 NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
3628 {
3629  if (_preconditioner.get() != nullptr)
3630  mooseError("More than one active Preconditioner detected");
3631 
3632  _preconditioner = pc;
3633 }
3634 
3635 MoosePreconditioner const *
3637 {
3638  return _preconditioner.get();
3639 }
3640 
3641 void
3643 {
3644  _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
3645 }
3646 
3647 void
3649  const std::set<MooseVariable *> & damped_vars)
3650 {
3651  for (const auto & var : damped_vars)
3652  var->computeIncrementAtQps(*_increment_vec);
3653 }
3654 
3655 void
3657  const std::set<MooseVariable *> & damped_vars)
3658 {
3659  for (const auto & var : damped_vars)
3660  var->computeIncrementAtNode(*_increment_vec);
3661 }
3662 
3663 void
3664 NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
3665 {
3666  // Obtain all blocks and variables covered by all kernels
3667  std::set<SubdomainID> input_subdomains;
3668  std::set<std::string> kernel_variables;
3669 
3670  bool global_kernels_exist = false;
3671  global_kernels_exist |= _scalar_kernels.hasActiveObjects();
3672  global_kernels_exist |= _nodal_kernels.hasActiveObjects();
3673 
3674  _kernels.subdomainsCovered(input_subdomains, kernel_variables);
3675  _dg_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3676  _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3677  _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3678  _constraints.subdomainsCovered(input_subdomains, kernel_variables);
3679 
3680 #ifdef MOOSE_KOKKOS_ENABLED
3681  _kokkos_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3682  _kokkos_nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
3683 #endif
3684 
3685  if (_fe_problem.haveFV())
3686  {
3687  std::vector<FVElementalKernel *> fv_elemental_kernels;
3689  .query()
3690  .template condition<AttribSystem>("FVElementalKernel")
3691  .queryInto(fv_elemental_kernels);
3692 
3693  for (auto fv_kernel : fv_elemental_kernels)
3694  {
3695  if (fv_kernel->blockRestricted())
3696  for (auto block_id : fv_kernel->blockIDs())
3697  input_subdomains.insert(block_id);
3698  else
3699  global_kernels_exist = true;
3700  kernel_variables.insert(fv_kernel->variable().name());
3701 
3702  // Check for lagrange multiplier
3703  if (dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel))
3704  kernel_variables.insert(dynamic_cast<FVScalarLagrangeMultiplierConstraint *>(fv_kernel)
3705  ->lambdaVariable()
3706  .name());
3707  }
3708 
3709  std::vector<FVFluxKernel *> fv_flux_kernels;
3711  .query()
3712  .template condition<AttribSystem>("FVFluxKernel")
3713  .queryInto(fv_flux_kernels);
3714 
3715  for (auto fv_kernel : fv_flux_kernels)
3716  {
3717  if (fv_kernel->blockRestricted())
3718  for (auto block_id : fv_kernel->blockIDs())
3719  input_subdomains.insert(block_id);
3720  else
3721  global_kernels_exist = true;
3722  kernel_variables.insert(fv_kernel->variable().name());
3723  }
3724 
3725  std::vector<FVInterfaceKernel *> fv_interface_kernels;
3727  .query()
3728  .template condition<AttribSystem>("FVInterfaceKernel")
3729  .queryInto(fv_interface_kernels);
3730 
3731  for (auto fvik : fv_interface_kernels)
3732  if (auto scalar_fvik = dynamic_cast<FVScalarLagrangeMultiplierInterface *>(fvik))
3733  kernel_variables.insert(scalar_fvik->lambdaVariable().name());
3734 
3735  std::vector<FVFluxBC *> fv_flux_bcs;
3737  .query()
3738  .template condition<AttribSystem>("FVFluxBC")
3739  .queryInto(fv_flux_bcs);
3740 
3741  for (auto fvbc : fv_flux_bcs)
3742  if (auto scalar_fvbc = dynamic_cast<FVBoundaryScalarLagrangeMultiplierConstraint *>(fvbc))
3743  kernel_variables.insert(scalar_fvbc->lambdaVariable().name());
3744  }
3745 
3746  // Check kernel coverage of subdomains (blocks) in your mesh
3747  if (!global_kernels_exist)
3748  {
3749  std::set<SubdomainID> difference;
3750  std::set_difference(mesh_subdomains.begin(),
3751  mesh_subdomains.end(),
3752  input_subdomains.begin(),
3753  input_subdomains.end(),
3754  std::inserter(difference, difference.end()));
3755 
3756  // there supposed to be no kernels on this lower-dimensional subdomain
3757  for (const auto & id : _mesh.interiorLowerDBlocks())
3758  difference.erase(id);
3759  for (const auto & id : _mesh.boundaryLowerDBlocks())
3760  difference.erase(id);
3761 
3762  if (!difference.empty())
3763  {
3764  std::vector<SubdomainID> difference_vec =
3765  std::vector<SubdomainID>(difference.begin(), difference.end());
3766  std::vector<SubdomainName> difference_names = _mesh.getSubdomainNames(difference_vec);
3767  std::stringstream missing_block_names;
3768  std::copy(difference_names.begin(),
3769  difference_names.end(),
3770  std::ostream_iterator<std::string>(missing_block_names, " "));
3771  std::stringstream missing_block_ids;
3772  std::copy(difference.begin(),
3773  difference.end(),
3774  std::ostream_iterator<unsigned int>(missing_block_ids, " "));
3775 
3776  mooseError("Each subdomain must contain at least one Kernel.\nThe following block(s) lack an "
3777  "active kernel: " +
3778  missing_block_names.str(),
3779  " (ids: ",
3780  missing_block_ids.str(),
3781  ")");
3782  }
3783  }
3784 
3785  // Check kernel use of variables
3786  std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
3787 
3788  std::set<VariableName> difference;
3789  std::set_difference(variables.begin(),
3790  variables.end(),
3791  kernel_variables.begin(),
3792  kernel_variables.end(),
3793  std::inserter(difference, difference.end()));
3794 
3795  // skip checks for varaibles defined on lower-dimensional subdomain
3796  std::set<VariableName> vars(difference);
3797  for (auto & var_name : vars)
3798  {
3799  auto blks = getSubdomainsForVar(var_name);
3800  for (const auto & id : blks)
3801  if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
3802  difference.erase(var_name);
3803  }
3804 
3805  if (!difference.empty())
3806  {
3807  std::stringstream missing_kernel_vars;
3808  std::copy(difference.begin(),
3809  difference.end(),
3810  std::ostream_iterator<std::string>(missing_kernel_vars, " "));
3811  mooseError("Each variable must be referenced by at least one active Kernel.\nThe following "
3812  "variable(s) lack an active kernel: " +
3813  missing_kernel_vars.str());
3814  }
3815 }
3816 
3817 bool
3819 {
3820  auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3821 
3822  return time_kernels.hasActiveObjects();
3823 }
3824 
3825 std::vector<std::string>
3827 {
3828  std::vector<std::string> variable_names;
3829  const auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
3830  if (time_kernels.hasActiveObjects())
3831  for (const auto & kernel : time_kernels.getObjects())
3832  variable_names.push_back(kernel->variable().name());
3833 
3834  return variable_names;
3835 }
3836 
3837 bool
3839 {
3840  // IntegratedBCs are for now the only objects we consider to be consuming
3841  // matprops on boundaries.
3842  if (_integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid))
3843  for (const auto & bc : _integrated_bcs.getActiveBoundaryObjects(bnd_id, tid))
3844  if (std::static_pointer_cast<MaterialPropertyInterface>(bc)->getMaterialPropertyCalled())
3845  return true;
3846 
3847  // Thin layer heat transfer in the heat_transfer module is being used on a boundary even though
3848  // it's an interface kernel. That boundary is external, on both sides of a gap in a mesh
3850  for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3851  if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3852  return true;
3853 
3854  // Because MortarConstraints do not inherit from BoundaryRestrictable, they are not sorted
3855  // by boundary in the MooseObjectWarehouse. So for now, we return true for all boundaries
3856  // Note: constraints are not threaded at this time
3857  if (_constraints.hasActiveObjects(/*tid*/ 0))
3858  for (const auto & ct : _constraints.getActiveObjects(/*tid*/ 0))
3859  if (auto mpi = std::dynamic_pointer_cast<MaterialPropertyInterface>(ct);
3860  mpi && mpi->getMaterialPropertyCalled())
3861  return true;
3862  return false;
3863 }
3864 
3865 bool
3867 {
3868  // InterfaceKernels are for now the only objects we consider to be consuming matprops on internal
3869  // boundaries.
3871  for (const auto & ik : _interface_kernels.getActiveBoundaryObjects(bnd_id, tid))
3872  if (std::static_pointer_cast<MaterialPropertyInterface>(ik)->getMaterialPropertyCalled())
3873  return true;
3874  return false;
3875 }
3876 
3877 bool
3879 {
3880  // DGKernels are for now the only objects we consider to be consuming matprops on
3881  // internal sides.
3882  if (_dg_kernels.hasActiveBlockObjects(subdomain_id, tid))
3883  for (const auto & dg : _dg_kernels.getActiveBlockObjects(subdomain_id, tid))
3884  if (std::static_pointer_cast<MaterialPropertyInterface>(dg)->getMaterialPropertyCalled())
3885  return true;
3886  // NOTE:
3887  // HDG kernels do not require face material properties on internal sides at this time.
3888  // The idea is to have element locality of HDG for hybridization
3889  return false;
3890 }
3891 
3892 bool
3894 {
3895  return _doing_dg;
3896 }
3897 
3898 void
3900 {
3903 }
3904 
3905 void
3907  const std::set<TagID> & vector_tags,
3908  const std::set<TagID> & matrix_tags)
3909 {
3910  parallel_object_only();
3911 
3912  try
3913  {
3914  for (auto & map_pr : _undisplaced_mortar_functors)
3915  map_pr.second(compute_type, vector_tags, matrix_tags);
3916 
3917  for (auto & map_pr : _displaced_mortar_functors)
3918  map_pr.second(compute_type, vector_tags, matrix_tags);
3919  }
3920  catch (MetaPhysicL::LogicError &)
3921  {
3922  mooseError(
3923  "We caught a MetaPhysicL error in NonlinearSystemBase::mortarConstraints. This is very "
3924  "likely due to AD not having a sufficiently large derivative container size. Please run "
3925  "MOOSE configure with the '--with-derivative-size=<n>' option");
3926  }
3927 }
3928 
3929 void
3931 {
3932  if (_auto_scaling_initd)
3933  return;
3934 
3935  // Want the libMesh count of variables, not MOOSE, e.g. I don't care about array variable counts
3936  const auto n_vars = system().n_vars();
3937 
3938  if (_scaling_group_variables.empty())
3939  {
3940  _var_to_group_var.reserve(n_vars);
3942 
3943  for (const auto var_number : make_range(n_vars))
3944  _var_to_group_var.emplace(var_number, var_number);
3945  }
3946  else
3947  {
3948  std::set<unsigned int> var_numbers, var_numbers_covered, var_numbers_not_covered;
3949  for (const auto var_number : make_range(n_vars))
3950  var_numbers.insert(var_number);
3951 
3953 
3954  for (const auto group_index : index_range(_scaling_group_variables))
3955  for (const auto & var_name : _scaling_group_variables[group_index])
3956  {
3957  if (!hasVariable(var_name) && !hasScalarVariable(var_name))
3958  mooseError("'",
3959  var_name,
3960  "', provided to the 'scaling_group_variables' parameter, does not exist in "
3961  "the nonlinear system.");
3962 
3963  const MooseVariableBase & var =
3964  hasVariable(var_name)
3965  ? static_cast<MooseVariableBase &>(getVariable(0, var_name))
3966  : static_cast<MooseVariableBase &>(getScalarVariable(0, var_name));
3967  auto map_pair = _var_to_group_var.emplace(var.number(), group_index);
3968  if (!map_pair.second)
3969  mooseError("Variable ", var_name, " is contained in multiple scaling grouplings");
3970  var_numbers_covered.insert(var.number());
3971  }
3972 
3973  std::set_difference(var_numbers.begin(),
3974  var_numbers.end(),
3975  var_numbers_covered.begin(),
3976  var_numbers_covered.end(),
3977  std::inserter(var_numbers_not_covered, var_numbers_not_covered.begin()));
3978 
3979  _num_scaling_groups = _scaling_group_variables.size() + var_numbers_not_covered.size();
3980 
3981  auto index = static_cast<unsigned int>(_scaling_group_variables.size());
3982  for (auto var_number : var_numbers_not_covered)
3983  _var_to_group_var.emplace(var_number, index++);
3984  }
3985 
3986  _variable_autoscaled.resize(n_vars, true);
3987  const auto & number_to_var_map = _vars[0].numberToVariableMap();
3988 
3990  for (const auto i : index_range(_variable_autoscaled))
3993  libmesh_map_find(number_to_var_map, i)->name()) !=
3995  _variable_autoscaled[i] = false;
3996 
3997  _auto_scaling_initd = true;
3998 }
3999 
4000 bool
4002 {
4004  return true;
4005 
4006  _console << "\nPerforming automatic scaling calculation\n" << std::endl;
4007 
4008  TIME_SECTION("computeScaling", 3, "Computing Automatic Scaling");
4009 
4010  // It's funny but we need to assemble our vector of scaling factors here otherwise we will be
4011  // applying scaling factors of 0 during Assembly of our scaling Jacobian
4013 
4014  // container for repeated access of element global dof indices
4015  std::vector<dof_id_type> dof_indices;
4016 
4017  if (!_auto_scaling_initd)
4018  setupScalingData();
4019 
4020  std::vector<Real> inverse_scaling_factors(_num_scaling_groups, 0);
4021  std::vector<Real> resid_inverse_scaling_factors(_num_scaling_groups, 0);
4022  std::vector<Real> jac_inverse_scaling_factors(_num_scaling_groups, 0);
4023  auto & dof_map = dofMap();
4024 
4025  // what types of scaling do we want?
4026  bool jac_scaling = _resid_vs_jac_scaling_param < 1. - TOLERANCE;
4027  bool resid_scaling = _resid_vs_jac_scaling_param > TOLERANCE;
4028 
4029  const NumericVector<Number> & scaling_residual = RHS();
4030 
4031  if (jac_scaling)
4032  {
4033  // if (!_auto_scaling_initd)
4034  // We need to reinit this when the number of dofs changes
4035  // but there is no good way to track that
4036  // In theory, it is the job of libmesh system to track this,
4037  // but this special matrix is not owned by libMesh system
4038  // Let us reinit eveytime since it is not expensive
4039  {
4040  auto init_vector = NumericVector<Number>::build(this->comm());
4041  init_vector->init(system().n_dofs(), system().n_local_dofs(), /*fast=*/false, PARALLEL);
4042 
4043  _scaling_matrix->clear();
4044  _scaling_matrix->init(*init_vector);
4045  }
4046 
4048  // Dispatch to derived classes to ensure that we use the correct matrix tag
4051  }
4052 
4053  if (resid_scaling)
4054  {
4057  // Dispatch to derived classes to ensure that we use the correct vector tag
4061  }
4062 
4063  // Did something bad happen during residual/Jacobian scaling computation?
4065  return false;
4066 
4067  auto examine_dof_indices = [this,
4068  jac_scaling,
4069  resid_scaling,
4070  &dof_map,
4071  &jac_inverse_scaling_factors,
4072  &resid_inverse_scaling_factors,
4073  &scaling_residual](const auto & dof_indices, const auto var_number)
4074  {
4075  for (auto dof_index : dof_indices)
4076  if (dof_map.local_index(dof_index))
4077  {
4078  if (jac_scaling)
4079  {
4080  // For now we will use the diagonal for determining scaling
4081  auto mat_value = (*_scaling_matrix)(dof_index, dof_index);
4082  auto & factor = jac_inverse_scaling_factors[_var_to_group_var[var_number]];
4083  factor = std::max(factor, std::abs(mat_value));
4084  }
4085  if (resid_scaling)
4086  {
4087  auto vec_value = scaling_residual(dof_index);
4088  auto & factor = resid_inverse_scaling_factors[_var_to_group_var[var_number]];
4089  factor = std::max(factor, std::abs(vec_value));
4090  }
4091  }
4092  };
4093 
4094  // Compute our scaling factors for the spatial field variables
4095  for (const auto & elem : _fe_problem.getCurrentAlgebraicElementRange())
4096  for (const auto i : make_range(system().n_vars()))
4098  {
4099  dof_map.dof_indices(elem, dof_indices, i);
4100  examine_dof_indices(dof_indices, i);
4101  }
4102 
4103  for (const auto i : make_range(system().n_vars()))
4104  if (_variable_autoscaled[i] && system().variable_type(i).family == SCALAR)
4105  {
4106  dof_map.SCALAR_dof_indices(dof_indices, i);
4107  examine_dof_indices(dof_indices, i);
4108  }
4109 
4110  if (resid_scaling)
4111  _communicator.max(resid_inverse_scaling_factors);
4112  if (jac_scaling)
4113  _communicator.max(jac_inverse_scaling_factors);
4114 
4115  if (jac_scaling && resid_scaling)
4116  for (MooseIndex(inverse_scaling_factors) i = 0; i < inverse_scaling_factors.size(); ++i)
4117  {
4118  // Be careful not to take log(0)
4119  if (!resid_inverse_scaling_factors[i])
4120  {
4121  if (!jac_inverse_scaling_factors[i])
4122  inverse_scaling_factors[i] = 1;
4123  else
4124  inverse_scaling_factors[i] = jac_inverse_scaling_factors[i];
4125  }
4126  else if (!jac_inverse_scaling_factors[i])
4127  // We know the resid is not zero
4128  inverse_scaling_factors[i] = resid_inverse_scaling_factors[i];
4129  else
4130  inverse_scaling_factors[i] =
4131  std::exp(_resid_vs_jac_scaling_param * std::log(resid_inverse_scaling_factors[i]) +
4132  (1 - _resid_vs_jac_scaling_param) * std::log(jac_inverse_scaling_factors[i]));
4133  }
4134  else if (jac_scaling)
4135  inverse_scaling_factors = jac_inverse_scaling_factors;
4136  else if (resid_scaling)
4137  inverse_scaling_factors = resid_inverse_scaling_factors;
4138  else
4139  mooseError("We shouldn't be calling this routine if we're not performing any scaling");
4140 
4141  // We have to make sure that our scaling values are not zero
4142  for (auto & scaling_factor : inverse_scaling_factors)
4143  if (scaling_factor == 0)
4144  scaling_factor = 1;
4145 
4146  // Now flatten the group scaling factors to the individual variable scaling factors
4147  std::vector<Real> flattened_inverse_scaling_factors(system().n_vars());
4148  for (const auto i : index_range(flattened_inverse_scaling_factors))
4149  flattened_inverse_scaling_factors[i] = inverse_scaling_factors[_var_to_group_var[i]];
4150 
4151  // Now set the scaling factors for the variables
4152  applyScalingFactors(flattened_inverse_scaling_factors);
4154  displaced_problem->systemBaseNonlinear(number()).applyScalingFactors(
4155  flattened_inverse_scaling_factors);
4156 
4157  _computed_scaling = true;
4158  return true;
4159 }
4160 
4161 void
4163 {
4164  if (!hasVector("scaling_factors"))
4165  // No variables have indicated they need scaling
4166  return;
4167 
4168  auto & scaling_vector = getVector("scaling_factors");
4169 
4170  const auto & lm_mesh = _mesh.getMesh();
4171  const auto & dof_map = dofMap();
4172 
4173  const auto & field_variables = _vars[0].fieldVariables();
4174  const auto & scalar_variables = _vars[0].scalars();
4175 
4176  std::vector<dof_id_type> dof_indices;
4177 
4178  for (const Elem * const elem :
4179  as_range(lm_mesh.active_local_elements_begin(), lm_mesh.active_local_elements_end()))
4180  for (const auto * const field_var : field_variables)
4181  {
4182  const auto & factors = field_var->arrayScalingFactor();
4183  for (const auto i : make_range(field_var->count()))
4184  {
4185  dof_map.dof_indices(elem, dof_indices, field_var->number() + i);
4186  for (const auto dof : dof_indices)
4187  scaling_vector.set(dof, factors[i]);
4188  }
4189  }
4190 
4191  for (const auto * const scalar_var : scalar_variables)
4192  {
4193  mooseAssert(scalar_var->count() == 1,
4194  "Scalar variables should always have only one component.");
4195  dof_map.SCALAR_dof_indices(dof_indices, scalar_var->number());
4196  for (const auto dof : dof_indices)
4197  scaling_vector.set(dof, scalar_var->scalingFactor());
4198  }
4199 
4200  // Parallel assemble
4201  scaling_vector.close();
4202 
4203  if (auto * displaced_problem = _fe_problem.getDisplacedProblem().get())
4204  // copy into the corresponding displaced system vector because they should be the exact same
4205  displaced_problem->systemBaseNonlinear(number()).getVector("scaling_factors") = scaling_vector;
4206 }
4207 
4208 bool
4210 {
4211  // Clear the iteration counters
4212  _current_l_its.clear();
4213  _current_nl_its = 0;
4214 
4215  // Initialize the solution vector using a predictor and known values from nodal bcs
4217 
4218  // Now that the initial solution has ben set, potentially perform a residual/Jacobian evaluation
4219  // to determine variable scaling factors
4220  if (_automatic_scaling)
4221  {
4222  const bool scaling_succeeded = computeScaling();
4223  if (!scaling_succeeded)
4224  return false;
4225  }
4226 
4227  // We do not know a priori what variable a global degree of freedom corresponds to, so we need a
4228  // map from global dof to scaling factor. We just use a ghosted NumericVector for that mapping
4230 
4231  return true;
4232 }
4233 
4234 void
4236 {
4237  if (matrixFromColoring())
4238  LibmeshPetscCall(MatFDColoringDestroy(&_fdcoloring));
4239 }
4240 
4243 {
4244  if (!_fsp)
4245  mooseError("No field split preconditioner is present for this system");
4246 
4247  return *_fsp;
4248 }
std::string name(const ElemQuality q)
std::vector< std::shared_ptr< TimeIntegrator > > _time_integrators
Time integrator.
Definition: SystemBase.h:1049
virtual void setSolutionUDotDotOld(const NumericVector< Number > &u_dotdot_old)
virtual void reinitNeighborPhys(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points, const THREAD_ID tid)=0
virtual void residualSetup(THREAD_ID tid=0) const
NumericVector< Number > & getResidualTimeVector()
Return a numeric vector that is associated with the time tag.
void setActiveMaterialProperties(const std::unordered_set< unsigned int > &mat_prop_ids, const THREAD_ID tid)
Record and set the material properties required by the current computing thread.
MooseObjectTagWarehouse< NodalKernelBase > _nodal_kernels
NodalKernels for each thread.
NumericVector< Number > * _Re_time
residual vector for time contributions
void computeJacobianBlocks(std::vector< JacobianBlock *> &blocks)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:42
TagID _Re_time_tag
Tag for time contribution residual.
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
unsigned int size(THREAD_ID tid=0) const
Return how many kernels we store in the current warehouse.
std::map< std::pair< BoundaryID, BoundaryID >, PenetrationLocator * > _penetration_locators
dof_id_type end_dof(const processor_id_type proc) const
bool empty() const
virtual void addKernel(const std::string &kernel_name, const std::string &name, InputParameters &parameters)
Adds a kernel.
virtual void setSolutionUDotDot(const NumericVector< Number > &udotdot)
Set transient term used by residual and Jacobian evaluation.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
A kernel for hybridized finite element formulations.
Definition: HDGKernel.h:17
MooseObjectTagWarehouse< ResidualObject > _kokkos_nodal_kernels
KOKKOS_INLINE_FUNCTION const T * find(const T &target, const T *const begin, const T *const end)
Find a value in an array.
Definition: KokkosUtils.h:30
void reinitIncrementAtNodeForDampers(THREAD_ID tid, const std::set< MooseVariable *> &damped_vars)
Compute the incremental change in variables at nodes for dampers.
void overwriteNodeFace(NumericVector< Number > &soln)
Called from explicit time stepping to overwrite boundary positions (explicit dynamics).
Base class for deriving general dampers.
Definition: GeneralDamper.h:21
bool _use_pre_smo_residual
Whether to use the pre-SMO initial residual in the relative convergence check.
MoosePreconditioner const * getPreconditioner() const
virtual const char * what() const
Get out the error message.
std::vector< std::pair< MooseVariableFEBase *, MooseVariableFEBase * > > & couplingEntries(const THREAD_ID tid, const unsigned int nl_sys_num)
void findImplicitGeometricCouplingEntries(GeometricSearchData &geom_search_data, std::unordered_map< dof_id_type, std::vector< dof_id_type >> &graph)
Finds the implicit sparsity graph between geometrically related dofs.
void setupDampers()
Setup damping stuff (called before we actually start)
Real _initial_residual
The initial (i.e., 0th nonlinear iteration) residual, see setPreSMOResidual for a detailed explanatio...
unsigned int n_threads()
std::vector< bool > _variable_autoscaled
Container to hold flag if variable is to participate in autoscaling.
void zeroVectorForResidual(const std::string &vector_name)
virtual void cacheResidualNeighbor(const THREAD_ID tid) override
bool identifyVariableGroupsInNL() const
Whether to identify variable groups in nonlinear systems.
void zeroTaggedVectors(const std::set< TagID > &tags)
Zero all vectors for given tags.
Definition: SystemBase.C:693
const std::set< SubdomainID > & interiorLowerDBlocks() const
Definition: MooseMesh.h:1429
Base class for split-based preconditioners.
Definition: Split.h:25
bool hasActiveBlockObjects(THREAD_ID tid=0) const
std::shared_ptr< DisplacedProblem > displaced_problem
bool hasVector(const std::string &tag_name) const
Check if the named vector exists in the system.
Definition: SystemBase.C:924
virtual void checkExceptionAndStopSolve(bool print_message=true)
Check to see if an exception has occurred on any processor and, if possible, force the solve to fail...
unsigned int n_comp(const unsigned int s, const unsigned int var) const
void applyScalingFactors(const std::vector< Real > &inverse_scaling_factors)
Applies scaling factors to the system&#39;s variables.
Definition: SystemBase.C:1495
void computeNodalBCsResidualAndJacobian()
compute the residual and Jacobian for nodal boundary conditions
bool _debugging_residuals
true if debugging residuals
SCALAR
BoundaryID _secondary_boundary
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
const std::map< SubdomainID, std::vector< std::shared_ptr< T > > > & getActiveBlockObjects(THREAD_ID tid=0) const
MooseObjectTagWarehouse< ResidualObject > _kokkos_kernels
unsigned int TagID
Definition: MooseTypes.h:210
Real computeDamping(const NumericVector< Number > &solution, const NumericVector< Number > &update)
Compute damping.
virtual void reinitNode(const Node *node, const THREAD_ID tid) override
virtual void setPreviousNewtonSolution(const NumericVector< Number > &soln)
virtual void predictorCleanup(NumericVector< libMesh::Number > &ghosted_solution)
Perform cleanup tasks after application of predictor to solution vector.
bool _assemble_constraints_separately
Whether or not to assemble the residual and Jacobian after the application of each constraint...
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:3153
auto exp(const T &)
TagID systemMatrixTag() const override
Return the Matrix Tag ID for System.
MPI_Info info
NumericVector< Number > & solution()
Definition: SystemBase.h:196
bool hasObjects(THREAD_ID tid=0) const
Convenience functions for determining if objects exist.
MooseObjectTagWarehouse< DGKernelBase > _dg_kernels
virtual bool haveFV() const override
returns true if this problem includes/needs finite volume functionality.
char ** blocks
face_info_iterator ownedFaceInfoBegin()
Iterators to owned faceInfo objects.
Definition: MooseMesh.C:1548
PARALLEL
void reinitIncrementAtQpsForDampers(THREAD_ID tid, const std::set< MooseVariable *> &damped_vars)
Compute the incremental change in variables at QPs for dampers.
void computeResidualAndJacobianInternal(const std::set< TagID > &vector_tags, const std::set< TagID > &matrix_tags)
Compute residual and Jacobian from contributions not related to constraints, such as nodal boundary c...
void addImplicitGeometricCouplingEntriesToJacobian(bool add=true)
If called with true this will add entries into the jacobian to link together degrees of freedom that ...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
bool areCoupled(const unsigned int ivar, const unsigned int jvar, const unsigned int nl_sys_num) const
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
unsigned int number() const
Get variable number coming from libMesh.
virtual void initialSetup() override
Setup Functions.
Data structure used to hold penetration information.
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Combine two vector parameters into a single vector of pairs.
const std::vector< std::shared_ptr< NodalConstraint > > & getActiveNodalConstraints() const
Access methods for active objects.
NumericVector< Number > * _u_dot_old
old solution vector for u^dot
Definition: SystemBase.h:1011
bool _has_nodalbc_diag_save_in
If there is a nodal BC having diag_save_in.
virtual void reinitScalars(const THREAD_ID tid, bool reinit_for_derivative_reordering=false) override
fills the VariableValue arrays for scalar variables from the solution vector
Base class for automatic differentiation Dirichlet BCs.
void add(std::shared_ptr< MooseObject > obj)
add adds a new object to the warehouse and stores attributes/metadata about it for running queries/fi...
Definition: TheWarehouse.C:116
virtual void getDiracElements(std::set< const Elem *> &elems) override
Fills "elems" with the elements that should be looped over for Dirac Kernels.
void setupDM()
Setup the PETSc DM object (when appropriate)
void setCurrentlyComputingResidual(bool currently_computing_residual) final
Set whether or not the problem is in the process of computing the residual.
void checkKernelCoverage(const std::set< SubdomainID > &mesh_subdomains) const
void addDGKernel(std::string dg_kernel_name, const std::string &name, InputParameters &parameters)
Adds a DG kernel.
void computeJacobian(libMesh::SparseMatrix< Number > &jacobian, const std::set< TagID > &tags)
Associate jacobian to systemMatrixTag, and then form a matrix for all the tags.
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseBase.h:131
char ** vars
virtual void reinitNeighborPhys(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points, const THREAD_ID tid) override
virtual TagID addVectorTag(const TagName &tag_name, const Moose::VectorTagType type=Moose::VECTOR_TAG_RESIDUAL)
Create a Tag.
Definition: SubProblem.C:92
std::vector< T * > & queryInto(std::vector< T *> &results, Args &&... args)
queryInto executes the query and stores the results in the given vector.
Definition: TheWarehouse.h:311
void getNodeDofs(dof_id_type node_id, std::vector< dof_id_type > &dofs)
Base class for all Constraint types.
Definition: Constraint.h:19
std::set< TagID > _nl_vector_tags
Vector tags to temporarily store all tags associated with the current system.
std::vector< std::string > _ignore_variables_for_autoscaling
A container for variables that do not partipate in autoscaling.
void residualSetup() override
Base boundary condition of a Dirichlet type.
virtual void associateVectorToTag(NumericVector< Number > &vec, TagID tag)
Associate a vector for a given tag.
Definition: SystemBase.C:981
virtual void reinitNodes(const std::vector< dof_id_type > &nodes, const THREAD_ID tid) override
MooseObjectTagWarehouse< NodalBCBase > _nodal_bcs
virtual void addJacobianOffDiagScalar(unsigned int ivar, const THREAD_ID tid=0)
virtual void setException(const std::string &message)
Set an exception, which is stored at this point by toggling a member variable in this class...
const std::vector< std::shared_ptr< NodeFaceConstraint > > & getActiveNodeFaceConstraints(BoundaryID boundary_id, bool displaced) const
virtual void customSetup(const ExecFlagType &exec_type, THREAD_ID tid=0) const
Factory & _factory
Definition: SystemBase.h:989
MooseObjectWarehouse< T > & getVectorTagsObjectWarehouse(const std::set< TagID > &tags, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object at least has one of the given vector ta...
virtual const Node * queryNodePtr(const dof_id_type i) const
Definition: MooseMesh.C:875
const ElementPairInfo & getElemPairInfo(std::pair< const Elem *, const Elem *> elem_pair) const
MooseObjectWarehouseBase< Split > _splits
Decomposition splits.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
NumericVector< Number > * _u_dotdot
solution vector for u^dotdot
Definition: SystemBase.h:1008
bool hasActiveMortarConstraints(const std::pair< BoundaryID, BoundaryID > &mortar_interface_key, bool displaced) const
const Parallel::Communicator & comm() const
NumericVector< Number > & add_vector(std::string_view vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Solving a linear problem.
Definition: MooseTypes.h:849
bool _has_nodalbc_save_in
If there is a nodal BC having save_in.
bool hasActiveNodalConstraints() const
Deterimine if active objects exist.
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
Definition: SystemBase.h:360
MooseObjectWarehouse< NodalDamper > _nodal_dampers
Nodal Dampers for each thread.
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
virtual libMesh::NonlinearSolver< Number > * nonlinearSolver()=0
dof_id_type n_dofs(const unsigned int vn) const
Real preSMOResidual() const
The pre-SMO residual.
std::unique_ptr< libMesh::DiagonalMatrix< Number > > _scaling_matrix
A diagonal matrix used for computing scaling.
virtual void setSolutionUDotOld(const NumericVector< Number > &u_dot_old)
virtual void associateMatrixToTag(libMesh::SparseMatrix< Number > &matrix, TagID tag)
Associate a matrix to a tag.
Definition: SystemBase.C:1076
bool hasDiagSaveIn() const
Weather or not the nonlinear system has diagonal Jacobian save-ins.
This class provides an interface for common operations on field variables of both FE and FV types wit...
const Parallel::Communicator & _communicator
void updateActive(THREAD_ID tid=0) override
Update the various active lists.
Real initialResidual() const
The initial residual.
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
bool needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
Indicated whether this system needs material properties on interfaces.
const libMesh::ConstElemRange & getCurrentAlgebraicElementRange()
These are the element and nodes that contribute to the jacobian and residual for this local processor...
MooseObjectWarehouse< ResidualObject > _kokkos_preset_nodal_bcs
std::size_t _num_scaling_groups
The number of scaling groups.
const libMesh::ConstNodeRange & getCurrentAlgebraicNodeRange()
void computingScalingJacobian(bool computing_scaling_jacobian)
Setter for whether we&#39;re computing the scaling jacobian.
dof_id_type n_local_dofs(const unsigned int vn) const
virtual void updateActive(THREAD_ID tid=0) override
Update the active status of Kernels.
bool has_dofs(const unsigned int s=libMesh::invalid_uint) const
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
Data structure of nodes and their associated penetration information.
bool hasActiveNodeElemConstraints(SubdomainID secondary_id, SubdomainID primary_id, bool displaced) const
Real _pre_smo_residual
The pre-SMO residual, see setPreSMOResidual for a detailed explanation.
bool hasDampers()
Whether or not this system has dampers.
const std::vector< std::shared_ptr< NodeElemConstraintBase > > & getActiveNodeElemConstraints(SubdomainID secondary_id, SubdomainID primary_id, bool displaced) const
bool _compute_scaling_once
Whether the scaling factors should only be computed once at the beginning of the simulation through a...
void computeResidualTags(const std::set< TagID > &tags)
Form multiple tag-associated residual vectors for all the given tags.
virtual void cacheJacobianNeighbor(const THREAD_ID tid) override
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
bool _has_save_in
If there is any Kernel or IntegratedBC having save_in.
virtual const Node & nodeRef(const dof_id_type i) const
Definition: MooseMesh.C:849
TagID _Ke_system_tag
Tag for system contribution Jacobian.
virtual void setResidual(NumericVector< libMesh::Number > &residual, const THREAD_ID tid) override
dof_id_type n_dofs() const
virtual void disassociateMatrixFromTag(libMesh::SparseMatrix< Number > &matrix, TagID tag)
Disassociate a matrix from a tag.
Definition: SystemBase.C:1088
Scope guard for starting and stopping Floating Point Exception Trapping.
virtual void zero()=0
auto max(const L &left, const R &right)
Serves as a base class for DGKernel and ADDGKernel.
Definition: DGKernelBase.h:32
const Variable & variable(const unsigned int c) const override
void constraintResiduals(NumericVector< Number > &residual, bool displaced)
Add residual contributions from Constraints.
Base class for MOOSE preconditioners.
NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const libMesh::ParallelType type)
Adds a solution length vector to the system.
virtual GeometricSearchData & geomSearchData() override
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor > _undisplaced_mortar_functors
Functors for computing undisplaced mortar constraints.
Specialization for filling multiple "small" preconditioning matrices simulatenously.
virtual bool matrixFromColoring() const
Whether a system matrix is formed from coloring.
Definition: SolverSystem.h:102
void update()
Update the system (doing libMesh magic)
Definition: SystemBase.C:1243
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override
void addScalarKernel(const std::string &kernel_name, const std::string &name, InputParameters &parameters)
Adds a scalar kernel.
void addBoundaryCondition(const std::string &bc_name, const std::string &name, InputParameters &parameters)
Adds a boundary condition.
bool hasActiveNodeFaceConstraints(BoundaryID boundary_id, bool displaced) const
void computeResidual(NumericVector< Number > &residual, TagID tag_id)
Form a residual vector for a given tag.
virtual void addCachedResidualDirectly(NumericVector< libMesh::Number > &residual, const THREAD_ID tid)
Allows for all the residual contributions that are currently cached to be added directly into the vec...
virtual unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:891
EXTERN_C_BEGIN PetscErrorCode DMCreate_Moose(DM)
bool hasActiveBoundaryObjects(THREAD_ID tid=0) const
bool _need_residual_ghosted
Whether or not a ghosted copy of the residual needs to be made.
virtual void activateAllMatrixTags()
Make all existing matrices active.
Definition: SystemBase.C:1131
virtual const std::string & name() const
Definition: SystemBase.C:1340
virtual void jacobianSetup()
Definition: SystemBase.C:1606
virtual bool containsTimeKernel() override
If the system has a kernel that corresponds to a time derivative.
void onTimestepBegin()
Called at the beginning of the time step.
const ConstBndNodeRange & getCurrentAlgebraicBndNodeRange()
std::vector< dof_id_type, Threads::scalable_allocator< dof_id_type > > Row
MooseObjectTagWarehouse< DiracKernelBase > _dirac_kernels
Dirac Kernel storage for each thread.
void closeTaggedMatrices(const std::set< TagID > &tags)
Close all matrices associated the tags.
Definition: SystemBase.C:1060
TagID _Re_non_time_tag
Tag for non-time contribution residual.
std::set< TagID > _nl_matrix_tags
Matrix tags to temporarily store all tags associated with the current system.
void solutionInvalidAccumulation()
Pass the number of solution invalid occurrences from current iteration to cumulative counters...
void set_basic_system_only()
std::map< std::pair< BoundaryID, BoundaryID >, NearestNodeLocator * > _nearest_node_locators
void setPredictor(std::shared_ptr< Predictor > predictor)
Real _resid_vs_jac_scaling_param
The param that indicates the weighting of the residual vs the Jacobian in determining variable scalin...
bool _auto_scaling_initd
Whether we&#39;ve initialized the automatic scaling data structures.
virtual void computeScalingResidual()=0
Compute a "residual" for automatic scaling purposes.
bool _doing_dg
true if DG is active (optimization reasons)
virtual libMesh::DofMap & dofMap()
Gets writeable reference to the dof map.
Definition: SystemBase.C:1163
void syncIteration()
Sync iteration counts to main processor.
virtual void deactivateAllMatrixTags()
Make matrices inactive.
Definition: SystemBase.C:1119
unsigned int number() const
void computeResidualAndJacobianTags(const std::set< TagID > &vector_tags, const std::set< TagID > &matrix_tags)
Form possibly multiple tag-associated vectors and matrices.
SERIAL
bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
Indicated whether this system needs material properties on boundaries.
MooseObjectWarehouse< DirichletBCBase > _preset_nodal_bcs
std::unordered_map< unsigned int, unsigned int > _var_to_group_var
A map from variable index to group variable index and it&#39;s associated (inverse) scaling factor...
const std::vector< std::shared_ptr< T > > & getActiveObjects(THREAD_ID tid=0) const
Retrieve complete vector to the active all/block/boundary restricted objects for a given thread...
std::vector< unsigned int > _current_l_its
virtual void zero()=0
virtual std::unique_ptr< Base > create()=0
void setCurrentNonlinearSystem(const unsigned int nl_sys_num)
std::shared_ptr< T > getActiveObject(const std::string &name, THREAD_ID tid=0) const
This is the common base class for the three main kernel types implemented in MOOSE, Kernel, VectorKernel and ArrayKernel.
Definition: KernelBase.h:23
dof_id_type id() const
std::vector< dof_id_type > _secondary_nodes
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3488
void min(const T &r, T &o, Request &req) const
void computeDiracContributions(const std::set< TagID > &tags, bool is_jacobian)
void updateActive(THREAD_ID tid)
Update active objects of Warehouses owned by NonlinearSystemBase.
TheWarehouse & theWarehouse() const
unsigned int n_vars
void reinitMaterialsNeighbor(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful=true, const std::deque< MaterialBase *> *reinit_mats=nullptr)
reinit materials on the neighboring element face
const ElementPairList & getElemPairs() const
FieldSplitPreconditionerBase & getFieldSplitPreconditioner()
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor > _displaced_mortar_functors
Functors for computing displaced mortar constraints.
GHOSTED
boundary_id_type BoundaryID
Real referenceResidual() const
The reference residual used in relative convergence check.
void addSplit(const std::string &split_name, const std::string &name, InputParameters &parameters)
Adds a split.
bool _automatic_scaling
Whether to automatically scale the variables.
Definition: SystemBase.h:1055
void computeJacobianTags(const std::set< TagID > &tags)
Computes multiple (tag associated) Jacobian matricese.
SolutionInvalidity & solutionInvalidity()
Get the SolutionInvalidity for this app.
Definition: MooseApp.h:179
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
virtual void timestepSetup(THREAD_ID tid=0) const
NonlinearSystemBase & currentNonlinearSystem()
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
void computingScalingResidual(bool computing_scaling_residual)
Setter for whether we&#39;re computing the scaling residual.
virtual void setupDM()=0
setup the data management data structure that manages the field split
SubProblem & subproblem()
Definition: SystemBase.h:101
std::vector< std::string > _vecs_to_zero_for_residual
vectors that will be zeroed before a residual computation
virtual TagID addMatrixTag(TagName tag_name)
Create a Tag.
Definition: SubProblem.C:311
std::shared_ptr< Split > getSplit(const std::string &name)
Retrieves a split by name.
std::unique_ptr< NumericVector< Number > > solution
SubProblem & _subproblem
The subproblem for whom this class holds variable data, etc; this can either be the governing finite ...
Definition: SystemBase.h:983
void subdomainsCovered(std::set< SubdomainID > &subdomains_covered, std::set< std::string > &unique_variables, THREAD_ID tid=0) const
Populates a set of covered subdomains and the associated variable names.
void addImplicitGeometricCouplingEntries(GeometricSearchData &geom_search_data)
Adds entries to the Jacobian in the correct positions for couplings coming from dofs being coupled th...
virtual void addHDGKernel(const std::string &kernel_name, const std::string &name, InputParameters &parameters)
Adds a hybridized discontinuous Galerkin (HDG) kernel.
virtual void zero_rows(std::vector< numeric_index_type > &rows, T diag_value=0.0)
MooseObjectTagWarehouse< KernelBase > _kernels
MooseObjectWarehouse< T > & getMatrixTagObjectWarehouse(TagID tag_id, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has the given matrix tag...
bool shouldEvaluatePreSMOResidual() const
We offer the option to check convergence against the pre-SMO residual.
std::vector< VectorTag > getVectorTags(const std::set< TagID > &tag_ids) const
Definition: SubProblem.C:172
Moose::SolveType _type
Definition: SolverParams.h:19
Base class for deriving nodal dampers.
Definition: NodalDamper.h:27
virtual void disassociateVectorFromTag(NumericVector< Number > &vec, TagID tag)
Disassociate a given vector from a given tag.
virtual void cacheResidual(const THREAD_ID tid) override
void mooseDeprecated(Args &&... args)
Emit a deprecated code/feature message with the given stringified, concatenated args.
Definition: MooseError.h:374
virtual void computeScalingJacobian()=0
Compute a "Jacobian" for automatic scaling purposes.
bool errorOnJacobianNonzeroReallocation() const
Will return True if the user wants to get an error when a nonzero is reallocated in the Jacobian by P...
This is the ElementPairLocator class.
This is the ElementPairInfo class.
std::map< BoundaryID, std::shared_ptr< ElementPairLocator > > _element_pair_locators
Moose::CouplingType coupling() const
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1157
void setupScalingData()
Setup group scaling containers.
virtual GeometricSearchData & geomSearchData()=0
const bool & usePreSMOResidual() const
Whether we are using pre-SMO residual in relative convergence checks.
MooseObjectTagWarehouse< HDGKernel > _hybridized_kernels
const std::set< SubdomainID > & boundaryLowerDBlocks() const
Definition: MooseMesh.h:1433
std::vector< SubdomainName > getSubdomainNames(const std::vector< SubdomainID > &subdomain_ids) const
Get the associated subdomainNames for the subdomain ids that are passed in.
Definition: MooseMesh.C:1807
auto log(const T &)
AuxiliarySystem & getAuxiliarySystem()
virtual void initialSetup(THREAD_ID tid=0) const
Convenience methods for calling object setup methods.
virtual void updateGeomSearch(GeometricSearchData::GeometricSearchType type=GeometricSearchData::ALL) override
void closeTaggedVectors(const std::set< TagID > &tags)
Close all vectors for given tags.
Definition: SystemBase.C:667
Base class for deriving any boundary condition that works at nodes.
Definition: NodalBCBase.h:26
void computeResidualInternal(const std::set< TagID > &tags)
Compute the residual for a given tag.
virtual void prepareAssembly(const THREAD_ID tid) override
bool computeScaling()
Method used to obtain scaling factors for variables.
Interface for objects interacting with the PerfGraph.
virtual std::map< TagName, TagID > & getMatrixTags()
Return all matrix tags in the system, where a tag is represented by a map from name to ID...
Definition: SubProblem.h:253
virtual bool hasVariable(const std::string &var_name) const
Query a system for a variable.
Definition: SystemBase.C:851
virtual void setCurrentSubdomainID(const Elem *elem, const THREAD_ID tid) override
virtual void clearDiracInfo() override
Gets called before Dirac Kernels are asked to add the points they are supposed to be evaluated in...
void destroyColoring()
Destroy the coloring object if it exists.
virtual void reinitNodesNeighbor(const std::vector< dof_id_type > &nodes, const THREAD_ID tid) override
bool identify_variable_groups() const
virtual void close()=0
virtual void jacobianSetup(THREAD_ID tid=0) const
ConstraintWarehouse _constraints
Constraints storage object.
virtual void turnOffJacobian()
Turn off the Jacobian (must be called before equation system initialization)
void computeKokkosJacobian(const std::set< TagID > &tags)
Compute Jacobian with Kokkos objects.
TagID residualVectorTag() const override
const std::map< BoundaryID, std::vector< std::shared_ptr< T > > > & getActiveBoundaryObjects(THREAD_ID tid=0) const
void computingNonlinearResid(bool computing_nonlinear_residual) final
Set whether or not the problem is in the process of computing the nonlinear residual.
virtual void customSetup(const ExecFlagType &exec_type) override
ComputeType
The type of nonlinear computation being performed.
Definition: MooseTypes.h:781
Base class for deriving element dampers.
Definition: ElementDamper.h:33
Base interface for field split preconditioner.
bool _add_implicit_geometric_coupling_entries_to_jacobian
Whether or not to add implicit geometric couplings to the Jacobian for FDP.
const_iterator end() const
MooseObjectTagWarehouse< ResidualObject > _kokkos_integrated_bcs
virtual void addNodalKernel(const std::string &kernel_name, const std::string &name, InputParameters &parameters)
Adds a NodalKernel.
virtual unsigned int numMatrixTags() const
The total number of tags.
Definition: SubProblem.h:248
Base class for creating new types of boundary conditions.
tbb::split split
MooseApp & _app
Definition: SystemBase.h:988
FEProblemBase & _fe_problem
the governing finite element/volume problem
Definition: SystemBase.h:986
virtual void reinitElemPhys(const Elem *elem, const std::vector< Point > &phys_points_in_elem, const THREAD_ID tid)=0
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a scalar variable with specified number.
Definition: SystemBase.C:145
ParallelType type() const
virtual NumericVector< Number > & RHS()=0
std::vector< VariableWarehouse > _vars
Variable warehouses (one for each thread)
Definition: SystemBase.h:996
Provides a way for users to bail out of the current solve.
unsigned int _n_residual_evaluations
Total number of residual evaluations that have been performed.
virtual void update()
void addDiracKernel(const std::string &kernel_name, const std::string &name, InputParameters &parameters)
Adds a Dirac kernel.
virtual void close()=0
void computeJacobianInternal(const std::set< TagID > &tags)
Form multiple matrices for all the tags.
bool _has_diag_save_in
If there is any Kernel or IntegratedBC having diag_save_in.
Base class for creating new types of nodal kernels.
const FEType & variable_type(const unsigned int i) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual std::shared_ptr< const DisplacedProblem > getDisplacedProblem() const
std::unique_ptr< NumericVector< Number > > _residual_copy
Copy of the residual vector, or nullptr if a copy is not needed.
virtual void subdomainSetup()
Definition: SystemBase.C:1592
bool hasSaveIn() const
Weather or not the nonlinear system has save-ins.
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:78
TagID timeVectorTag() const override
Ideally, we should not need this API.
Class for containing MooseEnum item information.
Definition: MooseEnumItem.h:18
MooseMesh & _mesh
Definition: SystemBase.h:991
virtual std::vector< std::string > timeKernelVariableNames() override
Returns the names of the variables that have time derivative kernels in the system.
bool hasActiveObjects(THREAD_ID tid=0) const
NonlinearSystemBase(FEProblemBase &problem, libMesh::System &sys, const std::string &name)
void max(const T &r, T &o, Request &req) const
NumericVector< Number > * _u_dot
solution vector for u^dot
Definition: SystemBase.h:1006
MooseObjectWarehouse< ElementDamper > _element_dampers
Element Dampers for each thread.
void addObject(std::shared_ptr< Constraint > object, THREAD_ID tid=0, bool recurse=true) override
Add Constraint object to the warehouse.
virtual void setSolutionUDot(const NumericVector< Number > &udot)
Set transient term used by residual and Jacobian evaluation.
virtual libMesh::SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:1024
virtual void augmentSparsity(libMesh::SparsityPattern::Graph &sparsity, std::vector< dof_id_type > &n_nz, std::vector< dof_id_type > &n_oz) override
Will modify the sparsity pattern to add logical geometric connections.
const std::vector< std::shared_ptr< MortarConstraintBase > > & getActiveMortarConstraints(const std::pair< BoundaryID, BoundaryID > &mortar_interface_key, bool displaced) const
MooseObjectWarehouse< T > & getMatrixTagsObjectWarehouse(const std::set< TagID > &tags, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has one of the given matrix tags...
void computeScalarKernelsJacobians(const std::set< TagID > &tags)
void setInitialResidual(Real r)
Record the initial residual (for later relative convergence check)
Query query()
query creates and returns an initialized a query object for querying objects from the warehouse...
Definition: TheWarehouse.h:466
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num)=0
const_iterator begin() const
bool needInternalNeighborSideMaterial(SubdomainID subdomain_id, THREAD_ID tid) const
Indicates whether this system needs material properties on internal sides.
void computeKokkosResidual(const std::set< TagID > &tags)
Compute residual with Kokkos objects.
Base class for deriving dampers.
Definition: Damper.h:24
bool getFailNextNonlinearConvergenceCheck() const
Whether it will skip further residual evaluations and fail the next nonlinear convergence check(s) ...
Base class shared by AD and non-AD scalar kernels.
void addDamper(const std::string &damper_name, const std::string &name, InputParameters &parameters)
Adds a damper.
IntRange< T > make_range(T beg, T end)
virtual MooseMesh & mesh() override
virtual void postAddResidualObject(ResidualObject &)
Called after any ResidualObject-derived objects are added to the system.
const std::vector< VariableName > & getVariableNames() const
Definition: SystemBase.h:860
virtual void preInit() override
This is called prior to the libMesh system has been init&#39;d.
Definition: SolverSystem.C:32
virtual void updateActive(THREAD_ID tid=0) override
Update the active status of Kernels.
virtual void timestepSetup() override
bool hasActiveElemElemConstraints(const InterfaceID interface_id, bool displaced) const
bool _off_diagonals_in_auto_scaling
Whether to include off diagonals when determining automatic scaling factors.
virtual unsigned int numVectorTags(const Moose::VectorTagType type=Moose::VECTOR_TAG_ANY) const
The total number of tags, which can be limited to the tag type.
Definition: SubProblem.C:195
void reinitNodeFace(const Node &secondary_node, const BoundaryID secondary_boundary, const PenetrationInfo &info, const bool displaced)
Reinitialize quantities such as variables, residuals, Jacobians, materials for node-face constraints...
Base class for deriving any boundary condition of a integrated type.
SolverParams & solverParams(unsigned int solver_sys_num=0)
Get the solver parameters.
NumericVector< Number > * _residual_ghosted
ghosted form of the residual
void setCachedJacobian(GlobalDataKey)
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
Definition: Assembly.C:4492
TagID _Re_tag
Used for the residual vector from PETSc.
virtual void customSetup(const ExecFlagType &exec_type)
Definition: SystemBase.C:1585
MooseObjectWarehouse< T > & getVectorTagObjectWarehouse(TagID tag_id, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has the given vector tag...
void computeNodalBCs(NumericVector< Number > &residual)
Enforces nodal boundary conditions.
const std::set< SubdomainID > & getSubdomainsForVar(unsigned int var_number) const
Definition: SystemBase.h:761
virtual void addObject(std::shared_ptr< T > object, THREAD_ID tid=0, bool recurse=true)
Adds an object to the storage structure.
bool ignoreZerosInJacobian() const
Will return true if zeros in the Jacobian are to be dropped from the sparsity pattern.
const ExecFlagType EXEC_PRE_KERNELS
Definition: Moose.C:56
libMesh::System & _sys
void mortarConstraints(Moose::ComputeType compute_type, const std::set< TagID > &vector_tags, const std::set< TagID > &matrix_tags)
Do mortar constraint residual/jacobian computations.
NumericVector< Number > * _increment_vec
increment vector
InterfaceKernelBase is the base class for all InterfaceKernel type classes.
QueryCache & condition(Args &&... args)
Adds a new condition to the query.
Definition: TheWarehouse.h:284
bool doingDG() const
Getter for _doing_dg.
void computeResidualTag(NumericVector< Number > &residual, TagID tag_id)
Computes residual for a given tag.
void addConstraint(const std::string &c_name, const std::string &name, InputParameters &parameters)
Adds a Constraint.
virtual TagName vectorTagName(const TagID tag) const
Retrieve the name associated with a TagID.
Definition: SubProblem.C:221
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
bool hasKokkosObjects() const
bool restoreOriginalNonzeroPattern() const
face_info_iterator ownedFaceInfoEnd()
Definition: MooseMesh.C:1557
void constraintJacobians(const SparseMatrix< Number > &jacobian_to_view, bool displaced)
Add jacobian contributions from Constraints.
virtual libMesh::System & system() override
Get the reference to the libMesh system.
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:90
std::vector< BoundaryID > getBoundaryIDs(const Elem *const elem, const unsigned short int side) const
Returns a vector of boundary IDs for the requested element on the requested side. ...
bool preSolve()
Perform some steps to get ready for the solver.
void full_sparsity_pattern_needed()
const std::unordered_map< std::pair< BoundaryID, BoundaryID >, AutomaticMortarGeneration > & getMortarInterfaces(bool on_displaced) const
bool _has_constraints
Whether or not this system has any Constraints.
dof_id_type first_dof(const processor_id_type proc) const
bool _computed_scaling
Flag used to indicate whether we have already computed the scaling Jacobian.
void remove_algebraic_ghosting_functor(GhostingFunctor &evaluable_functor)
MooseObjectTagWarehouse< InterfaceKernelBase > _interface_kernels
unsigned int n_vars() const
void resetSolutionInvalidCurrentIteration()
Reset the number of solution invalid occurrences back to zero.
NumericVector< Number > & residualVector(TagID tag)
Return a residual vector that is associated with the residual tag.
NumericVector< Number > & solutionOld()
Definition: SystemBase.h:197
NumericVector< Number > & getResidualNonTimeVector()
Return a numeric vector that is associated with the nontime tag.
virtual bool hasScalarVariable(const std::string &var_name) const
Definition: SystemBase.C:876
const TagName PREVIOUS_NL_SOLUTION_TAG
Definition: MooseTypes.C:28
processor_id_type processor_id() const
std::shared_ptr< Predictor > _predictor
If predictor is active, this is non-NULL.
void enforceNodalConstraintsResidual(NumericVector< Number > &residual)
Enforce nodal constraints.
void subdomainsCovered(std::set< SubdomainID > &subdomains_covered, std::set< std::string > &unique_variables, THREAD_ID tid=0) const
Update supplied subdomain and variable coverate containters.
MooseObjectWarehouse< GeneralDamper > _general_dampers
General Dampers.
const std::vector< std::shared_ptr< ElemElemConstraint > > & getActiveElemElemConstraints(InterfaceID interface_id, bool displaced) const
virtual void initialSetup()
Setup Functions.
Definition: SystemBase.C:1558
bool defaultGhosting()
Whether or not the user has requested default ghosting ot be on.
Definition: SubProblem.h:144
void addCachedJacobian(GlobalDataKey)
Adds the values that have been cached by calling cacheJacobian() and or cacheJacobianNeighbor() to th...
Definition: Assembly.C:3815
virtual void cacheJacobian(const THREAD_ID tid) override
auto min(const L &left, const R &right)
void jacobianSetup() override
void setKokkosInitialSolution()
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector by name.
Definition: SystemBase.C:933
MooseObjectTagWarehouse< IntegratedBCBase > _integrated_bcs
const DofMap & get_dof_map() const
virtual void reinitNodeFace(const Node *node, BoundaryID bnd_id, const THREAD_ID tid) override
virtual void residualSetup()
Definition: SystemBase.C:1599
virtual void reinitOffDiagScalars(const THREAD_ID tid) override
processor_id_type processor_id() const
virtual void addObject(std::shared_ptr< T > object, THREAD_ID tid=0, bool recurse=true) override
Adds an object to the storage structure.
virtual void setNeighborSubdomainID(const Elem *elem, unsigned int side, const THREAD_ID tid) override
virtual void addResidualScalar(const THREAD_ID tid=0)
virtual void subdomainSetup(THREAD_ID tid=0) const
virtual void jacobianSetup() override
virtual void addCachedResidual(const THREAD_ID tid) override
std::vector< std::vector< std::string > > _scaling_group_variables
A container of variable groupings that can be used in scaling calculations.
void addInterfaceKernel(std::string interface_kernel_name, const std::string &name, InputParameters &parameters)
Adds an interface kernel.
FieldSplitPreconditionerBase * _fsp
The field split preconditioner if this sytem is using one.
auto index_range(const T &sizable)
virtual NumericVector< Number > & residualCopy() override
void reinitMortarFunctors()
Update the mortar functors if the mesh has changed.
virtual NumericVector< Number > & residualGhosted() override
DiracKernelBase is the base class for all DiracKernel type classes.
virtual void updateActive(THREAD_ID tid=0)
Updates the active objects storage.
MooseObjectTagWarehouse< ResidualObject > _kokkos_nodal_bcs
void assembleScalingVector()
Assemble the numeric vector of scaling factors such that it can be used during assembly of the system...
BoundaryID _primary_boundary
NumericVector< Number > * _u_dotdot_old
old solution vector for u^dotdot
Definition: SystemBase.h:1013
Base variable class.
MooseObjectTagWarehouse< ScalarKernelBase > _scalar_kernels
virtual void residualEnd(THREAD_ID tid=0) const
void setConstraintSecondaryValues(NumericVector< Number > &solution, bool displaced)
Sets the value of constrained variables in the solution vector.
unsigned int THREAD_ID
Definition: MooseTypes.h:209
uint8_t dof_id_type
virtual void addJacobianScalar(const THREAD_ID tid=0)
NearestNodeLocator & _nearest_node
const std::map< dof_id_type, std::vector< dof_id_type > > & nodeToElemMap()
If not already created, creates a map from every node to all elements to which they are connected...
Definition: MooseMesh.C:1216
const std::set< SubdomainID > & meshSubdomains() const
Returns a read-only reference to the set of subdomains currently present in the Mesh.
Definition: MooseMesh.C:3211
ParallelType
virtual void addCachedJacobian(const THREAD_ID tid) override
virtual void timestepSetup()
Definition: SystemBase.C:1578
virtual ~NonlinearSystemBase()
virtual void preInit() override
This is called prior to the libMesh system has been init&#39;d.
virtual void residualSetup() override
void setPreconditioner(std::shared_ptr< MoosePreconditioner > pc)
Sets a preconditioner.
virtual void localize(std::vector< T > &v_local) const=0
MooseObjectWarehouse< ADDirichletBCBase > _ad_preset_nodal_bcs
virtual libMesh::System & system() override
Get the reference to the libMesh system.
Key structure for APIs manipulating global vectors/matrices.
Definition: Assembly.h:844