www.mooseframework.org
NonlinearSystemBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "NonlinearSystemBase.h"
11 #include "AuxiliarySystem.h"
12 #include "Problem.h"
13 #include "FEProblem.h"
14 #include "MooseVariableFE.h"
15 #include "MooseVariableScalar.h"
16 #include "PetscSupport.h"
17 #include "Factory.h"
18 #include "ParallelUniqueId.h"
19 #include "ThreadedElementLoop.h"
20 #include "MaterialData.h"
21 #include "ComputeResidualThread.h"
22 #include "ComputeJacobianThread.h"
25 #include "ComputeDiracThread.h"
32 #include "TimeKernel.h"
33 #include "BoundaryCondition.h"
34 #include "PresetNodalBC.h"
35 #include "NodalBCBase.h"
36 #include "IntegratedBCBase.h"
37 #include "DGKernel.h"
38 #include "InterfaceKernelBase.h"
39 #include "ElementDamper.h"
40 #include "NodalDamper.h"
41 #include "GeneralDamper.h"
42 #include "DisplacedProblem.h"
43 #include "NearestNodeLocator.h"
44 #include "PenetrationLocator.h"
45 #include "NodalConstraint.h"
46 #include "NodeFaceConstraint.h"
47 #include "NodeElemConstraint.h"
48 #include "MortarConstraint.h"
49 #include "ElemElemConstraint.h"
50 #include "ScalarKernel.h"
51 #include "Parser.h"
52 #include "Split.h"
54 #include "MooseMesh.h"
55 #include "MooseUtils.h"
56 #include "MooseApp.h"
57 #include "NodalKernel.h"
58 #include "DiracKernel.h"
59 #include "NodalKernel.h"
60 #include "TimeIntegrator.h"
61 #include "Predictor.h"
62 #include "Assembly.h"
63 #include "ElementPairLocator.h"
64 #include "ODETimeKernel.h"
67 #include "MaxVarNDofsPerElem.h"
68 #include "MaxVarNDofsPerNode.h"
69 #include "ADKernel.h"
70 #include "ADPresetNodalBC.h"
71 #include "Moose.h"
72 
73 // libMesh
74 #include "libmesh/nonlinear_solver.h"
75 #include "libmesh/quadrature_gauss.h"
76 #include "libmesh/dense_vector.h"
77 #include "libmesh/boundary_info.h"
78 #include "libmesh/petsc_matrix.h"
79 #include "libmesh/petsc_vector.h"
80 #include "libmesh/petsc_nonlinear_solver.h"
81 #include "libmesh/numeric_vector.h"
82 #include "libmesh/mesh.h"
83 #include "libmesh/dense_subvector.h"
84 #include "libmesh/dense_submatrix.h"
85 #include "libmesh/dof_map.h"
86 #include "libmesh/sparse_matrix.h"
87 #include "libmesh/petsc_matrix.h"
88 
89 // PETSc
90 #ifdef LIBMESH_HAVE_PETSC
91 #include "petscsnes.h"
92 #if !PETSC_VERSION_LESS_THAN(3, 3, 0)
93 #include <PetscDMMoose.h>
94 EXTERN_C_BEGIN
95 extern PetscErrorCode DMCreate_Moose(DM);
96 EXTERN_C_END
97 #endif
98 #endif
99 
101  System & sys,
102  const std::string & name)
103  : SystemBase(fe_problem, name, Moose::VAR_NONLINEAR),
104  ConsoleStreamInterface(fe_problem.getMooseApp()),
105  PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "NonlinearSystemBase"),
106  _fe_problem(fe_problem),
107  _sys(sys),
108  _last_rnorm(0.),
109  _last_nl_rnorm(0.),
110  _l_abs_step_tol(1e-10),
111  _initial_residual_before_preset_bcs(0.),
112  _initial_residual_after_preset_bcs(0.),
113  _current_nl_its(0),
114  _compute_initial_residual_before_preset_bcs(true),
115  _current_solution(NULL),
116  _residual_ghosted(NULL),
117  _serialized_solution(*NumericVector<Number>::build(_communicator).release()),
118  _solution_previous_nl(NULL),
119  _residual_copy(*NumericVector<Number>::build(_communicator).release()),
120  _u_dot(NULL),
121  _u_dotdot(NULL),
122  _u_dot_old(NULL),
123  _u_dotdot_old(NULL),
124  _Re_time_tag(-1),
125  _Re_time(NULL),
126  _Re_non_time_tag(-1),
127  _Re_non_time(NULL),
128  _scalar_kernels(/*threaded=*/false),
129  _nodal_bcs(/*threaded=*/false),
130  _preset_nodal_bcs(/*threaded=*/false),
131  _ad_preset_nodal_bcs(/*threaded=*/false),
132  _splits(/*threaded=*/false),
133  _increment_vec(NULL),
134  _pc_side(Moose::PCS_DEFAULT),
135  _ksp_norm(Moose::KSPN_UNPRECONDITIONED),
136  _use_finite_differenced_preconditioner(false),
137  _have_decomposition(false),
138  _use_field_split_preconditioner(false),
139  _add_implicit_geometric_coupling_entries_to_jacobian(false),
140  _assemble_constraints_separately(false),
141  _need_serialized_solution(false),
142  _need_residual_copy(false),
143  _need_residual_ghosted(false),
144  _debugging_residuals(false),
145  _doing_dg(false),
146  _n_iters(0),
147  _n_linear_iters(0),
148  _n_residual_evaluations(0),
149  _final_residual(0.),
150  _computing_initial_residual(false),
151  _print_all_var_norms(false),
152  _has_save_in(false),
153  _has_diag_save_in(false),
154  _has_nodalbc_save_in(false),
155  _has_nodalbc_diag_save_in(false),
156  _compute_residual_tags_timer(registerTimedSection("computeResidualTags", 5)),
157  _compute_residual_internal_timer(registerTimedSection("computeResidualInternal", 3)),
158  _kernels_timer(registerTimedSection("Kernels", 3)),
159  _scalar_kernels_timer(registerTimedSection("ScalarKernels", 3)),
160  _nodal_kernels_timer(registerTimedSection("NodalKernels", 3)),
161  _nodal_kernel_bcs_timer(registerTimedSection("NodalKernelBCs", 3)),
162  _nodal_bcs_timer(registerTimedSection("NodalBCs", 3)),
163  _compute_jacobian_tags_timer(registerTimedSection("computeJacobianTags", 5)),
164  _compute_jacobian_blocks_timer(registerTimedSection("computeJacobianBlocks", 3)),
165  _compute_dampers_timer(registerTimedSection("computeDampers", 3)),
166  _compute_dirac_timer(registerTimedSection("computeDirac", 3))
167 {
169  // Don't need to add the matrix - it already exists (for now)
171 
172  // The time matrix tag is not normally used - but must be added to the system
173  // in case it is so that objects can have 'time' in their matrix tags by default
174  _fe_problem.addMatrixTag("TIME");
175 
176  _Re_tag = _fe_problem.addVectorTag("RESIDUAL");
177 }
178 
180 {
181  delete &_serialized_solution;
182  delete &_residual_copy;
183 }
184 
185 void
187 {
188  if (_fe_problem.hasDampers())
189  setupDampers();
190 
191  _current_solution = _sys.current_local_solution.get();
192 
194  _serialized_solution.init(_sys.n_dofs(), false, SERIAL);
195 
197  _residual_copy.init(_sys.n_dofs(), false, SERIAL);
198 
199  Moose::perf_log.push("maxVarNDofsPerElem()", "Setup");
200  MaxVarNDofsPerElem mvndpe(_fe_problem, *this);
201  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mvndpe);
202  _max_var_n_dofs_per_elem = mvndpe.max();
203  _communicator.max(_max_var_n_dofs_per_elem);
204  auto displaced_problem = _fe_problem.getDisplacedProblem();
205  if (displaced_problem)
206  displaced_problem->nlSys().assignMaxVarNDofsPerElem(_max_var_n_dofs_per_elem);
207  Moose::perf_log.pop("maxVarNDofsPerElem()", "Setup");
208 
209  Moose::perf_log.push("maxVarNDofsPerNode()", "Setup");
210  MaxVarNDofsPerNode mvndpn(_fe_problem, *this);
211  Threads::parallel_reduce(*_mesh.getLocalNodeRange(), mvndpn);
212  _max_var_n_dofs_per_node = mvndpn.max();
213  _communicator.max(_max_var_n_dofs_per_node);
214  if (displaced_problem)
215  displaced_problem->nlSys().assignMaxVarNDofsPerNode(_max_var_n_dofs_per_node);
216  Moose::perf_log.pop("maxVarNDofsPerNode()", "Setup");
217 }
218 
219 void
221 {
223  _u_dot = &addVector("u_dot", true, GHOSTED);
225  _u_dot_old = &addVector("u_dot_old", true, GHOSTED);
227  _u_dotdot = &addVector("u_dotdot", true, GHOSTED);
229  _u_dotdot_old = &addVector("u_dotdot_old", true, GHOSTED);
230 }
231 
232 void
234 {
235  system().set_basic_system_only();
236  nonlinearSolver()->jacobian = NULL;
237 }
238 
239 void
241 {
243  _solution_previous_nl = &addVector("u_previous_newton", true, GHOSTED);
244 }
245 
246 void
248 {
249  // call parent
251  // and update _current_solution
252  _current_solution = _sys.current_local_solution.get();
253 }
254 
255 void
257 {
258  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
259  {
260  _kernels.initialSetup(tid);
264  if (_doing_dg)
270  }
275 
276  // go over mortar interfaces and construct functors
277  const auto & undisplaced_mortar_interfaces = _fe_problem.getMortarInterfaces(/*displaced=*/false);
278  for (const auto & mortar_interface : undisplaced_mortar_interfaces)
279  {
280  auto master_slave_boundary_pair = mortar_interface.first;
281  const auto & mortar_generation_object = mortar_interface.second;
282 
283  auto & mortar_constraints =
284  _constraints.getActiveMortarConstraints(master_slave_boundary_pair, /*displaced=*/false);
285 
287  master_slave_boundary_pair,
289  mortar_constraints, mortar_generation_object, _fe_problem));
291  master_slave_boundary_pair,
293  mortar_constraints, mortar_generation_object, _fe_problem));
294  }
295 
296  const auto & displaced_mortar_interfaces = _fe_problem.getMortarInterfaces(/*displaced=*/true);
297  for (const auto & mortar_interface : displaced_mortar_interfaces)
298  {
299  mooseAssert(_fe_problem.getDisplacedProblem(),
300  "Cannot create displaced mortar functors when the displaced problem is null");
301  auto master_slave_boundary_pair = mortar_interface.first;
302  const auto & mortar_generation_object = mortar_interface.second;
303 
304  auto & mortar_constraints =
305  _constraints.getActiveMortarConstraints(master_slave_boundary_pair, /*displaced=*/true);
306 
308  master_slave_boundary_pair,
310  mortar_constraints, mortar_generation_object, *_fe_problem.getDisplacedProblem()));
312  master_slave_boundary_pair,
314  mortar_constraints, mortar_generation_object, *_fe_problem.getDisplacedProblem()));
315  }
316 }
317 
318 void
320 {
321  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
322  {
323  _kernels.timestepSetup(tid);
327  if (_doing_dg)
333  }
338 }
339 
340 void
341 NonlinearSystemBase::setDecomposition(const std::vector<std::string> & splits)
342 {
344  if (splits.size() && splits.size() != 1)
345  mooseError("Only a single top-level split is allowed in a Problem's decomposition.");
346 
347  if (splits.size())
348  {
349  _decomposition_split = splits[0];
350  _have_decomposition = true;
351  }
352  else
353  _have_decomposition = false;
354 }
355 
356 void
358 {
359  if (!_have_decomposition)
360  return;
361 
362  std::shared_ptr<Split> top_split = getSplit(_decomposition_split);
363  top_split->setup();
364 }
365 
366 void
368  const std::string & name,
369  InputParameters parameters)
370 {
371  parameters.set<SystemBase *>("_sys") = this;
372 
373  std::shared_ptr<TimeIntegrator> ti = _factory.create<TimeIntegrator>(type, name, parameters);
374  _time_integrator = ti;
375 }
376 
377 void
378 NonlinearSystemBase::addKernel(const std::string & kernel_name,
379  const std::string & name,
380  InputParameters parameters)
381 {
382  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
383  {
384  // Create the kernel object via the factory and add to warehouse
385  std::shared_ptr<KernelBase> kernel =
386  _factory.create<KernelBase>(kernel_name, name, parameters, tid);
389  _ad_jacobian_kernels.addObject(kernel, tid);
390  else
391  _kernels.addObject(kernel, tid);
392  }
393 
394  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
395  _has_save_in = true;
396  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
397  _has_diag_save_in = true;
398 }
399 
400 void
401 NonlinearSystemBase::addNodalKernel(const std::string & kernel_name,
402  const std::string & name,
403  InputParameters parameters)
404 {
405  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
406  {
407  // Create the kernel object via the factory and add to the warehouse
408  std::shared_ptr<NodalKernel> kernel =
409  _factory.create<NodalKernel>(kernel_name, name, parameters, tid);
410  _nodal_kernels.addObject(kernel, tid);
411  }
412 
413  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
414  _has_save_in = true;
415  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
416  _has_diag_save_in = true;
417 }
418 
419 void
420 NonlinearSystemBase::addScalarKernel(const std::string & kernel_name,
421  const std::string & name,
422  InputParameters parameters)
423 {
424  std::shared_ptr<ScalarKernel> kernel =
425  _factory.create<ScalarKernel>(kernel_name, name, parameters);
426  _scalar_kernels.addObject(kernel);
427 }
428 
429 void
430 NonlinearSystemBase::addBoundaryCondition(const std::string & bc_name,
431  const std::string & name,
432  InputParameters parameters)
433 {
434  // ThreadID
435  THREAD_ID tid = 0;
436 
437  // Create the object
438  std::shared_ptr<BoundaryCondition> bc =
439  _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
440 
441  // Active BoundaryIDs for the object
442  const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
443  _vars[tid].addBoundaryVar(boundary_ids, &bc->variable());
444 
445  // Cast to the various types of BCs
446  std::shared_ptr<NodalBCBase> nbc = std::dynamic_pointer_cast<NodalBCBase>(bc);
447  std::shared_ptr<IntegratedBCBase> ibc = std::dynamic_pointer_cast<IntegratedBCBase>(bc);
448 
449  // NodalBCBase
450  if (nbc)
451  {
452  _nodal_bcs.addObject(nbc);
453  _vars[tid].addBoundaryVars(boundary_ids, nbc->getCoupledVars());
454 
455  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
456  _has_nodalbc_save_in = true;
457  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
459 
460  // PresetNodalBC
461  std::shared_ptr<PresetNodalBC> pnbc = std::dynamic_pointer_cast<PresetNodalBC>(bc);
462  if (pnbc)
464 
465  std::shared_ptr<ADPresetNodalBC<RESIDUAL>> adpnbc =
467  if (adpnbc)
469  }
470 
471  // IntegratedBCBase
472  else if (ibc)
473  {
474  _integrated_bcs.addObject(ibc, tid);
475  _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
476 
477  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
478  _has_save_in = true;
479  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
480  _has_diag_save_in = true;
481 
482  for (tid = 1; tid < libMesh::n_threads(); tid++)
483  {
484  // Create the object
485  bc = _factory.create<BoundaryCondition>(bc_name, name, parameters, tid);
486 
487  // Active BoundaryIDs for the object
488  const std::set<BoundaryID> & boundary_ids = bc->boundaryIDs();
489  _vars[tid].addBoundaryVar(boundary_ids, &bc->variable());
490 
491  ibc = std::static_pointer_cast<IntegratedBCBase>(bc);
492 
493  _integrated_bcs.addObject(ibc, tid);
494  _vars[tid].addBoundaryVars(boundary_ids, ibc->getCoupledVars());
495  }
496  }
497 
498  else
499  mooseError("Unknown BoundaryCondition type for object named ", bc->name());
500 }
501 
502 void
503 NonlinearSystemBase::addConstraint(const std::string & c_name,
504  const std::string & name,
505  InputParameters parameters)
506 {
507  std::shared_ptr<Constraint> constraint = _factory.create<Constraint>(c_name, name, parameters);
508  _constraints.addObject(constraint);
509 
510  if (constraint && constraint->addCouplingEntriesToJacobian())
512 }
513 
514 void
515 NonlinearSystemBase::addDiracKernel(const std::string & kernel_name,
516  const std::string & name,
517  InputParameters parameters)
518 {
519  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
520  {
521  std::shared_ptr<DiracKernel> kernel =
522  _factory.create<DiracKernel>(kernel_name, name, parameters, tid);
523  _dirac_kernels.addObject(kernel, tid);
524  }
525 }
526 
527 void
528 NonlinearSystemBase::addDGKernel(std::string dg_kernel_name,
529  const std::string & name,
530  InputParameters parameters)
531 {
532  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
533  {
534  auto dg_kernel = _factory.create<DGKernelBase>(dg_kernel_name, name, parameters, tid);
535  _dg_kernels.addObject(dg_kernel, tid);
536  }
537 
538  _doing_dg = true;
539 
540  if (parameters.get<std::vector<AuxVariableName>>("save_in").size() > 0)
541  _has_save_in = true;
542  if (parameters.get<std::vector<AuxVariableName>>("diag_save_in").size() > 0)
543  _has_diag_save_in = true;
544 }
545 
546 void
547 NonlinearSystemBase::addInterfaceKernel(std::string interface_kernel_name,
548  const std::string & name,
549  InputParameters parameters)
550 {
551  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
552  {
553  std::shared_ptr<InterfaceKernelBase> interface_kernel =
554  _factory.create<InterfaceKernelBase>(interface_kernel_name, name, parameters, tid);
555 
556  const std::set<BoundaryID> & boundary_ids = interface_kernel->boundaryIDs();
557  _vars[tid].addBoundaryVar(boundary_ids, &interface_kernel->variable());
558 
559  _interface_kernels.addObject(interface_kernel, tid);
560  _vars[tid].addBoundaryVars(boundary_ids, interface_kernel->getCoupledVars());
561  }
562 }
563 
564 void
565 NonlinearSystemBase::addDamper(const std::string & damper_name,
566  const std::string & name,
567  InputParameters parameters)
568 {
569  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
570  {
571  std::shared_ptr<Damper> damper = _factory.create<Damper>(damper_name, name, parameters, tid);
572 
573  // Attempt to cast to the damper types
574  std::shared_ptr<ElementDamper> ed = std::dynamic_pointer_cast<ElementDamper>(damper);
575  std::shared_ptr<NodalDamper> nd = std::dynamic_pointer_cast<NodalDamper>(damper);
576  std::shared_ptr<GeneralDamper> gd = std::dynamic_pointer_cast<GeneralDamper>(damper);
577 
578  if (gd)
579  {
581  break; // not threaded
582  }
583  else if (ed)
584  _element_dampers.addObject(ed, tid);
585  else if (nd)
586  _nodal_dampers.addObject(nd, tid);
587  else
588  mooseError("Invalid damper type");
589  }
590 }
591 
592 void
593 NonlinearSystemBase::addSplit(const std::string & split_name,
594  const std::string & name,
595  InputParameters parameters)
596 {
597  std::shared_ptr<Split> split = _factory.create<Split>(split_name, name, parameters);
599 }
600 
601 std::shared_ptr<Split>
602 NonlinearSystemBase::getSplit(const std::string & name)
603 {
604  return _splits.getActiveObject(name);
605 }
606 
607 void
608 NonlinearSystemBase::zeroVectorForResidual(const std::string & vector_name)
609 {
610  for (unsigned int i = 0; i < _vecs_to_zero_for_residual.size(); ++i)
611  if (vector_name == _vecs_to_zero_for_residual[i])
612  return;
613 
614  _vecs_to_zero_for_residual.push_back(vector_name);
615 }
616 
617 void
618 NonlinearSystemBase::computeResidualTag(NumericVector<Number> & residual, TagID tag_id)
619 {
620 
621  _nl_vector_tags.clear();
622  _nl_vector_tags.insert(tag_id);
624 
626 
628 
630 }
631 
632 void
633 NonlinearSystemBase::computeResidual(NumericVector<Number> & residual, TagID tag_id)
634 {
635  mooseDeprecated(" Please use computeResidualTag");
636 
637  computeResidualTag(residual, tag_id);
638 }
639 
640 void
641 NonlinearSystemBase::computeResidualTags(const std::set<TagID> & tags)
642 {
643  TIME_SECTION(_compute_residual_tags_timer);
644 
645  bool required_residual = tags.find(residualVectorTag()) == tags.end() ? false : true;
646 
648 
649  // not suppose to do anythin on matrix
651 
653 
654  for (const auto & numeric_vec : _vecs_to_zero_for_residual)
655  if (hasVector(numeric_vec))
656  {
657  NumericVector<Number> & vec = getVector(numeric_vec);
658  vec.close();
659  vec.zero();
660  }
661 
662  try
663  {
664  zeroTaggedVectors(tags);
666  closeTaggedVectors(tags);
667 
668  if (required_residual)
669  {
670  auto & residual = getVector(residualVectorTag());
671  if (_time_integrator)
672  _time_integrator->postResidual(residual);
673  else
674  residual += *_Re_non_time;
675  residual.close();
676  }
677  computeNodalBCs(tags);
678  closeTaggedVectors(tags);
679 
680  // If we are debugging residuals we need one more assignment to have the ghosted copy up to
681  // date
682  if (_need_residual_ghosted && _debugging_residuals && required_residual)
683  {
684  auto & residual = getVector(residualVectorTag());
685 
686  *_residual_ghosted = residual;
687  _residual_ghosted->close();
688  }
689  // Need to close and update the aux system in case residuals were saved to it.
692  if (hasSaveIn())
694  }
695  catch (MooseException & e)
696  {
697  // The buck stops here, we have already handled the exception by
698  // calling stopSolve(), it is now up to PETSc to return a
699  // "diverged" reason during the next solve.
700  }
701 
702  // not supposed to do anything on matrix
704 }
705 
706 void
708 {
709  if (_time_integrator)
710  _time_integrator->preSolve();
711  if (_predictor.get())
712  _predictor->timestepSetup();
713 }
714 
715 void
717 {
719 
720  NumericVector<Number> & initial_solution(solution());
721  if (_predictor.get() && _predictor->shouldApply())
722  {
723  _predictor->apply(initial_solution);
724  _fe_problem.predictorCleanup(initial_solution);
725  }
726 
727  // do nodal BC
729  for (const auto & bnode : bnd_nodes)
730  {
731  BoundaryID boundary_id = bnode->_bnd_id;
732  Node * node = bnode->_node;
733 
734  if (node->processor_id() == processor_id())
735  {
736  // reinit variables in nodes
737  _fe_problem.reinitNodeFace(node, boundary_id, 0);
738 
740  {
741  const auto & preset_bcs = _preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
742  for (const auto & preset_bc : preset_bcs)
743  preset_bc->computeValue(initial_solution);
744  }
746  {
747  const auto & preset_bcs_res = _ad_preset_nodal_bcs.getActiveBoundaryObjects(boundary_id);
748  for (const auto & preset_bc : preset_bcs_res)
749  preset_bc->computeValue(initial_solution);
750  }
751  }
752  }
753 
754  _sys.solution->close();
755  update();
756 
757  // Set constraint slave values
758  setConstraintSlaveValues(initial_solution, false);
759 
761  setConstraintSlaveValues(initial_solution, true);
762 }
763 
764 void
765 NonlinearSystemBase::setPredictor(std::shared_ptr<Predictor> predictor)
766 {
767  _predictor = predictor;
768 }
769 
770 void
772 {
773  _kernels.subdomainSetup(subdomain, tid);
774  _ad_jacobian_kernels.subdomainSetup(subdomain, tid);
775  _nodal_kernels.subdomainSetup(subdomain, tid);
776  _element_dampers.subdomainSetup(subdomain, tid);
777  _nodal_dampers.subdomainSetup(subdomain, tid);
778 }
779 
780 NumericVector<Number> &
782 {
783  if (!_Re_time)
784  {
786  _Re_time = &addVector(_Re_time_tag, false, GHOSTED);
787  }
788  return *_Re_time;
789 }
790 
791 NumericVector<Number> &
793 {
794  if (!_Re_non_time)
795  {
797  _Re_non_time = &addVector(_Re_non_time_tag, false, GHOSTED);
798  }
799  return *_Re_non_time;
800 }
801 
802 NumericVector<Number> &
804 {
805  mooseDeprecated("Please use getVector()");
806  switch (tag)
807  {
808  case 0:
809  return getResidualNonTimeVector();
810 
811  case 1:
812  return getResidualTimeVector();
813 
814  default:
815  mooseError("The required residual vector is not available");
816  }
817 }
818 
819 void
821 {
822  if (_time_integrator)
823  {
824  _time_integrator->preStep();
825  _time_integrator->computeTimeDerivatives();
826  }
827 }
828 
829 void
831 {
832  THREAD_ID tid = 0; // constraints are going to be done single-threaded
833  residual.close();
835  {
836  const auto & ncs = _constraints.getActiveNodalConstraints();
837  for (const auto & nc : ncs)
838  {
839  std::vector<dof_id_type> & slave_node_ids = nc->getSlaveNodeId();
840  std::vector<dof_id_type> & master_node_ids = nc->getMasterNodeId();
841 
842  if ((slave_node_ids.size() > 0) && (master_node_ids.size() > 0))
843  {
844  _fe_problem.reinitNodes(master_node_ids, tid);
845  _fe_problem.reinitNodesNeighbor(slave_node_ids, tid);
846  nc->computeResidual(residual);
847  }
848  }
849  _fe_problem.addCachedResidualDirectly(residual, tid);
850  residual.close();
851  }
852 }
853 
854 void
856 {
857  if (!hasMatrix(systemMatrixTag()))
858  mooseError(" A system matrix is required");
859 
860  auto & jacobian = getMatrix(systemMatrixTag());
861  THREAD_ID tid = 0; // constraints are going to be done single-threaded
862  jacobian.close();
864  {
865  const auto & ncs = _constraints.getActiveNodalConstraints();
866  for (const auto & nc : ncs)
867  {
868  std::vector<dof_id_type> & slave_node_ids = nc->getSlaveNodeId();
869  std::vector<dof_id_type> & master_node_ids = nc->getMasterNodeId();
870 
871  if ((slave_node_ids.size() > 0) && (master_node_ids.size() > 0))
872  {
873  _fe_problem.reinitNodes(master_node_ids, tid);
874  _fe_problem.reinitNodesNeighbor(slave_node_ids, tid);
875  nc->computeJacobian(jacobian);
876  }
877  }
879  jacobian.close();
880  }
881 }
882 
883 void
884 NonlinearSystemBase::setConstraintSlaveValues(NumericVector<Number> & solution, bool displaced)
885 {
886  std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators =
887  NULL;
888 
889  if (!displaced)
890  {
891  GeometricSearchData & geom_search_data = _fe_problem.geomSearchData();
892  penetration_locators = &geom_search_data._penetration_locators;
893  }
894  else
895  {
896  GeometricSearchData & displaced_geom_search_data =
897  _fe_problem.getDisplacedProblem()->geomSearchData();
898  penetration_locators = &displaced_geom_search_data._penetration_locators;
899  }
900 
901  bool constraints_applied = false;
902 
903  for (const auto & it : *penetration_locators)
904  {
905  PenetrationLocator & pen_loc = *(it.second);
906 
907  std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;
908 
909  BoundaryID slave_boundary = pen_loc._slave_boundary;
910 
911  if (_constraints.hasActiveNodeFaceConstraints(slave_boundary, displaced))
912  {
913  const auto & constraints =
914  _constraints.getActiveNodeFaceConstraints(slave_boundary, displaced);
915 
916  for (unsigned int i = 0; i < slave_nodes.size(); i++)
917  {
918  dof_id_type slave_node_num = slave_nodes[i];
919  Node & slave_node = _mesh.nodeRef(slave_node_num);
920 
921  if (slave_node.processor_id() == processor_id())
922  {
923  if (pen_loc._penetration_info[slave_node_num])
924  {
925  PenetrationInfo & info = *pen_loc._penetration_info[slave_node_num];
926 
927  const Elem * master_elem = info._elem;
928  unsigned int master_side = info._side_num;
929 
930  // reinit variables at the node
931  _fe_problem.reinitNodeFace(&slave_node, slave_boundary, 0);
932 
934 
935  std::vector<Point> points;
936  points.push_back(info._closest_point);
937 
938  // reinit variables on the master element's face at the contact point
939  _fe_problem.setNeighborSubdomainID(master_elem, 0);
940  _fe_problem.reinitNeighborPhys(master_elem, master_side, points, 0);
941 
942  for (const auto & nfc : constraints)
943  if (nfc->shouldApply())
944  {
945  constraints_applied = true;
946  nfc->computeSlaveValue(solution);
947  }
948  }
949  }
950  }
951  }
952  }
953 
954  // go over NodeELemConstraints
955  std::set<dof_id_type> unique_slave_node_ids;
956 
957  for (const auto & slave_id : _mesh.meshSubdomains())
958  {
959  for (const auto & master_id : _mesh.meshSubdomains())
960  {
961  if (_constraints.hasActiveNodeElemConstraints(slave_id, master_id, displaced))
962  {
963  const auto & constraints =
964  _constraints.getActiveNodeElemConstraints(slave_id, master_id, displaced);
965 
966  // get unique set of ids of all nodes on current block
967  unique_slave_node_ids.clear();
968  const MeshBase & meshhelper = _mesh.getMesh();
969  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(slave_id),
970  meshhelper.active_subdomain_elements_end(slave_id)))
971  {
972  for (auto & n : elem->node_ref_range())
973  unique_slave_node_ids.insert(n.id());
974  }
975 
976  for (auto slave_node_id : unique_slave_node_ids)
977  {
978  Node & slave_node = _mesh.nodeRef(slave_node_id);
979 
980  // check if slave node is on current processor
981  if (slave_node.processor_id() == processor_id())
982  {
983  // This reinits the variables that exist on the slave node
984  _fe_problem.reinitNodeFace(&slave_node, slave_id, 0);
985 
986  // This will set aside residual and jacobian space for the variables that have dofs
987  // on the slave node
989 
990  for (const auto & nec : constraints)
991  {
992  if (nec->shouldApply())
993  {
994  constraints_applied = true;
995  nec->computeSlaveValue(solution);
996  }
997  }
998  }
999  }
1000  }
1001  }
1002  }
1003 
1004  // See if constraints were applied anywhere
1005  _communicator.max(constraints_applied);
1006 
1007  if (constraints_applied)
1008  {
1009  solution.close();
1010  update();
1011  }
1012 }
1013 
1014 void
1015 NonlinearSystemBase::constraintResiduals(NumericVector<Number> & residual, bool displaced)
1016 {
1017  // Make sure the residual is in a good state
1018  residual.close();
1019 
1020  std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators =
1021  NULL;
1022 
1023  if (!displaced)
1024  {
1025  GeometricSearchData & geom_search_data = _fe_problem.geomSearchData();
1026  penetration_locators = &geom_search_data._penetration_locators;
1027  }
1028  else
1029  {
1030  GeometricSearchData & displaced_geom_search_data =
1031  _fe_problem.getDisplacedProblem()->geomSearchData();
1032  penetration_locators = &displaced_geom_search_data._penetration_locators;
1033  }
1034 
1035  bool constraints_applied;
1036  bool residual_has_inserted_values = false;
1038  constraints_applied = false;
1039  for (const auto & it : *penetration_locators)
1040  {
1042  {
1043  // Reset the constraint_applied flag before each new constraint, as they need to be
1044  // assembled separately
1045  constraints_applied = false;
1046  }
1047  PenetrationLocator & pen_loc = *(it.second);
1048 
1049  std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;
1050 
1051  BoundaryID slave_boundary = pen_loc._slave_boundary;
1052 
1053  if (_constraints.hasActiveNodeFaceConstraints(slave_boundary, displaced))
1054  {
1055  const auto & constraints =
1056  _constraints.getActiveNodeFaceConstraints(slave_boundary, displaced);
1057 
1058  for (unsigned int i = 0; i < slave_nodes.size(); i++)
1059  {
1060  dof_id_type slave_node_num = slave_nodes[i];
1061  Node & slave_node = _mesh.nodeRef(slave_node_num);
1062 
1063  if (slave_node.processor_id() == processor_id())
1064  {
1065  if (pen_loc._penetration_info[slave_node_num])
1066  {
1067  PenetrationInfo & info = *pen_loc._penetration_info[slave_node_num];
1068 
1069  const Elem * master_elem = info._elem;
1070  unsigned int master_side = info._side_num;
1071 
1072  // *These next steps MUST be done in this order!*
1073 
1074  // This reinits the variables that exist on the slave node
1075  _fe_problem.reinitNodeFace(&slave_node, slave_boundary, 0);
1076 
1077  // This will set aside residual and jacobian space for the variables that have dofs on
1078  // the slave node
1080 
1081  std::vector<Point> points;
1082  points.push_back(info._closest_point);
1083 
1084  // reinit variables on the master element's face at the contact point
1085  _fe_problem.setNeighborSubdomainID(master_elem, 0);
1086  _fe_problem.reinitNeighborPhys(master_elem, master_side, points, 0);
1087 
1088  for (const auto & nfc : constraints)
1089  if (nfc->shouldApply())
1090  {
1091  constraints_applied = true;
1092  nfc->computeResidual();
1093 
1094  if (nfc->overwriteSlaveResidual())
1095  {
1096  _fe_problem.setResidual(residual, 0);
1097  residual_has_inserted_values = true;
1098  }
1099  else
1102  }
1103  }
1104  }
1105  }
1106  }
1108  {
1109  // Make sure that slave contribution to master are assembled, and ghosts have been
1110  // exchanged, as current masters might become slaves on next iteration and will need to
1111  // contribute their former slaves' contributions to the future masters. See if constraints
1112  // were applied anywhere
1113  _communicator.max(constraints_applied);
1114 
1115  if (constraints_applied)
1116  {
1117  // If any of the above constraints inserted values in the residual, it needs to be
1118  // assembled before adding the cached residuals below.
1119  _communicator.max(residual_has_inserted_values);
1120  if (residual_has_inserted_values)
1121  {
1122  residual.close();
1123  residual_has_inserted_values = false;
1124  }
1126  residual.close();
1127 
1129  *_residual_ghosted = residual;
1130  }
1131  }
1132  }
1134  {
1135  _communicator.max(constraints_applied);
1136 
1137  if (constraints_applied)
1138  {
1139  // If any of the above constraints inserted values in the residual, it needs to be assembled
1140  // before adding the cached residuals below.
1141  _communicator.max(residual_has_inserted_values);
1142  if (residual_has_inserted_values)
1143  residual.close();
1144 
1146  residual.close();
1147 
1149  *_residual_ghosted = residual;
1150  }
1151  }
1152 
1153  mortarResidualConstraints(displaced);
1154 
1155  // go over element-element constraint interface
1156  std::map<unsigned int, std::shared_ptr<ElementPairLocator>> * element_pair_locators = nullptr;
1157 
1158  if (!displaced)
1159  {
1160  GeometricSearchData & geom_search_data = _fe_problem.geomSearchData();
1161  element_pair_locators = &geom_search_data._element_pair_locators;
1162  }
1163  else
1164  {
1165  GeometricSearchData & displaced_geom_search_data =
1166  _fe_problem.getDisplacedProblem()->geomSearchData();
1167  element_pair_locators = &displaced_geom_search_data._element_pair_locators;
1168  }
1169 
1170  THREAD_ID tid = 0;
1171  for (const auto & it : *element_pair_locators)
1172  {
1173  ElementPairLocator & elem_pair_loc = *(it.second);
1174 
1175  if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
1176  {
1177  // ElemElemConstraint objects
1178  const auto & _element_constraints =
1179  _constraints.getActiveElemElemConstraints(it.first, displaced);
1180 
1181  // go over pair elements
1182  const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
1183  elem_pair_loc.getElemPairs();
1184  for (const auto & pr : elem_pairs)
1185  {
1186  const Elem * elem1 = pr.first;
1187  const Elem * elem2 = pr.second;
1188 
1189  if (elem1->processor_id() != processor_id())
1190  continue;
1191 
1192  const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
1193 
1194  // for each element process constraints on the
1195  for (const auto & ec : _element_constraints)
1196  {
1197  _fe_problem.setCurrentSubdomainID(elem1, tid);
1199  _fe_problem.setNeighborSubdomainID(elem2, tid);
1201 
1202  ec->subProblem().prepareShapes(ec->variable().number(), tid);
1203  ec->subProblem().prepareNeighborShapes(ec->variable().number(), tid);
1204 
1205  ec->reinit(info);
1206  ec->computeResidual();
1209  }
1211  }
1212  }
1213  }
1214 
1215  // go over NodeELemConstraints
1216  std::set<dof_id_type> unique_slave_node_ids;
1217 
1218  constraints_applied = false;
1219  residual_has_inserted_values = false;
1220  for (const auto & slave_id : _mesh.meshSubdomains())
1221  {
1222  for (const auto & master_id : _mesh.meshSubdomains())
1223  {
1224  if (_constraints.hasActiveNodeElemConstraints(slave_id, master_id, displaced))
1225  {
1226  const auto & constraints =
1227  _constraints.getActiveNodeElemConstraints(slave_id, master_id, displaced);
1228 
1229  // get unique set of ids of all nodes on current block
1230  unique_slave_node_ids.clear();
1231  const MeshBase & meshhelper = _mesh.getMesh();
1232  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(slave_id),
1233  meshhelper.active_subdomain_elements_end(slave_id)))
1234  {
1235  for (auto & n : elem->node_ref_range())
1236  unique_slave_node_ids.insert(n.id());
1237  }
1238 
1239  for (auto slave_node_id : unique_slave_node_ids)
1240  {
1241  Node & slave_node = _mesh.nodeRef(slave_node_id);
1242  // check if slave node is on current processor
1243  if (slave_node.processor_id() == processor_id())
1244  {
1245  // This reinits the variables that exist on the slave node
1246  _fe_problem.reinitNodeFace(&slave_node, slave_id, 0);
1247 
1248  // This will set aside residual and jacobian space for the variables that have dofs
1249  // on the slave node
1251 
1252  for (const auto & nec : constraints)
1253  {
1254  if (nec->shouldApply())
1255  {
1256  constraints_applied = true;
1257  nec->computeResidual();
1258 
1259  if (nec->overwriteSlaveResidual())
1260  {
1261  _fe_problem.setResidual(residual, 0);
1262  residual_has_inserted_values = true;
1263  }
1264  else
1267  }
1268  }
1270  }
1271  }
1272  }
1273  }
1274  }
1275  _communicator.max(constraints_applied);
1276 
1277  if (constraints_applied)
1278  {
1279  // If any of the above constraints inserted values in the residual, it needs to be assembled
1280  // before adding the cached residuals below.
1281  _communicator.max(residual_has_inserted_values);
1282  if (residual_has_inserted_values)
1283  residual.close();
1284 
1286  residual.close();
1287 
1289  *_residual_ghosted = residual;
1290  }
1291 }
1292 
1293 void
1295 {
1296  TIME_SECTION(_compute_residual_internal_timer);
1297 
1298  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1299  {
1300  _kernels.residualSetup(tid);
1303  if (_doing_dg)
1309  }
1314 
1315  // reinit scalar variables
1316  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
1318 
1319  // residual contributions from the domain
1320  PARALLEL_TRY
1321  {
1322  TIME_SECTION(_kernels_timer);
1323 
1324  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
1325 
1327 
1328  Threads::parallel_reduce(elem_range, cr);
1329 
1330  unsigned int n_threads = libMesh::n_threads();
1331  for (unsigned int i = 0; i < n_threads;
1332  i++) // Add any cached residuals that might be hanging around
1334  }
1335  PARALLEL_CATCH;
1336 
1337  // residual contributions from the scalar kernels
1338  PARALLEL_TRY
1339  {
1340  // do scalar kernels (not sure how to thread this)
1342  {
1343  TIME_SECTION(_scalar_kernels_timer);
1344 
1345  MooseObjectWarehouse<ScalarKernel> * scalar_kernel_warehouse;
1346  // This code should be refactored once we can do tags for scalar
1347  // kernels
1348  // Should redo this based on Warehouse
1349  if (!tags.size() || tags.size() == _fe_problem.numVectorTags())
1350  scalar_kernel_warehouse = &_scalar_kernels;
1351  else if (tags.size() == 1)
1352  scalar_kernel_warehouse =
1353  &(_scalar_kernels.getVectorTagObjectWarehouse(*(tags.begin()), 0));
1354  else
1355  // scalar_kernels is not threading
1356  scalar_kernel_warehouse = &(_scalar_kernels.getVectorTagsObjectWarehouse(tags, 0));
1357 
1358  bool have_scalar_contributions = false;
1359  const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
1360  for (const auto & scalar_kernel : scalars)
1361  {
1362  scalar_kernel->reinit();
1363  const std::vector<dof_id_type> & dof_indices = scalar_kernel->variable().dofIndices();
1364  const DofMap & dof_map = scalar_kernel->variable().dofMap();
1365  const dof_id_type first_dof = dof_map.first_dof();
1366  const dof_id_type end_dof = dof_map.end_dof();
1367  for (dof_id_type dof : dof_indices)
1368  {
1369  if (dof >= first_dof && dof < end_dof)
1370  {
1371  scalar_kernel->computeResidual();
1372  have_scalar_contributions = true;
1373  break;
1374  }
1375  }
1376  }
1377  if (have_scalar_contributions)
1379  }
1380  }
1381  PARALLEL_CATCH;
1382 
1383  // residual contributions from Block NodalKernels
1384  PARALLEL_TRY
1385  {
1387  {
1388  TIME_SECTION(_nodal_kernels_timer);
1389 
1391 
1392  ConstNodeRange & range = *_mesh.getLocalNodeRange();
1393 
1394  if (range.begin() != range.end())
1395  {
1396  _fe_problem.reinitNode(*range.begin(), 0);
1397 
1398  Threads::parallel_reduce(range, cnk);
1399 
1400  unsigned int n_threads = libMesh::n_threads();
1401  for (unsigned int i = 0; i < n_threads;
1402  i++) // Add any cached residuals that might be hanging around
1404  }
1405  }
1406  }
1407  PARALLEL_CATCH;
1408 
1409  // residual contributions from boundary NodalKernels
1410  PARALLEL_TRY
1411  {
1413  {
1414  TIME_SECTION(_nodal_kernel_bcs_timer);
1415 
1417 
1418  ConstBndNodeRange & bnd_node_range = *_mesh.getBoundaryNodeRange();
1419 
1420  Threads::parallel_reduce(bnd_node_range, cnk);
1421 
1422  unsigned int n_threads = libMesh::n_threads();
1423  for (unsigned int i = 0; i < n_threads;
1424  i++) // Add any cached residuals that might be hanging around
1426  }
1427  }
1428  PARALLEL_CATCH;
1429 
1430  if (_need_residual_copy)
1431  {
1432  _Re_non_time->close();
1433  _Re_non_time->localize(_residual_copy);
1434  }
1435 
1437  {
1438  _Re_non_time->close();
1440  _residual_ghosted->close();
1441  }
1442 
1443  PARALLEL_TRY { computeDiracContributions(tags, false); }
1444  PARALLEL_CATCH;
1445 
1447  {
1448  PARALLEL_TRY { enforceNodalConstraintsResidual(*_Re_non_time); }
1449  PARALLEL_CATCH;
1450  _Re_non_time->close();
1451  }
1452 
1453  // Add in Residual contributions from Constraints
1455  {
1456  PARALLEL_TRY
1457  {
1458  // Undisplaced Constraints
1460 
1461  // Displaced Constraints
1464 
1467  }
1468  PARALLEL_CATCH;
1469  _Re_non_time->close();
1470  }
1471 }
1472 
1473 void
1474 NonlinearSystemBase::computeNodalBCs(NumericVector<Number> & residual)
1475 {
1476  _nl_vector_tags.clear();
1477 
1478  auto & tags = _fe_problem.getVectorTags();
1479  for (auto & tag : tags)
1480  _nl_vector_tags.insert(tag.second);
1481 
1483  computeNodalBCs(residual, _nl_vector_tags);
1485 }
1486 
1487 void
1488 NonlinearSystemBase::computeNodalBCs(NumericVector<Number> & residual, const std::set<TagID> & tags)
1489 {
1491 
1492  computeNodalBCs(tags);
1493 
1495 }
1496 
1497 void
1498 NonlinearSystemBase::computeNodalBCs(const std::set<TagID> & tags)
1499 {
1500  // We need to close the diag_save_in variables on the aux system before NodalBCBases clear the
1501  // dofs on boundary nodes
1502  if (_has_save_in)
1504 
1505  PARALLEL_TRY
1506  {
1507  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
1508 
1509  if (!bnd_nodes.empty())
1510  {
1511  TIME_SECTION(_nodal_bcs_timer);
1512 
1513  MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
1514 
1515  // Select nodal kernels
1516  if (tags.size() == _fe_problem.numVectorTags() || !tags.size())
1517  nbc_warehouse = &_nodal_bcs;
1518  else if (tags.size() == 1)
1519  nbc_warehouse = &(_nodal_bcs.getVectorTagObjectWarehouse(*(tags.begin()), 0));
1520  else
1521  nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0));
1522 
1523  for (const auto & bnode : bnd_nodes)
1524  {
1525  BoundaryID boundary_id = bnode->_bnd_id;
1526  Node * node = bnode->_node;
1527 
1528  if (node->processor_id() == processor_id())
1529  {
1530  // reinit variables in nodes
1531  _fe_problem.reinitNodeFace(node, boundary_id, 0);
1532  if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id))
1533  {
1534  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
1535  for (const auto & nbc : bcs)
1536  if (nbc->shouldApply())
1537  nbc->computeResidual();
1538  }
1539  }
1540  }
1541  }
1542  }
1543  PARALLEL_CATCH;
1544 
1545  if (_Re_time)
1546  _Re_time->close();
1547  _Re_non_time->close();
1548 }
1549 
1550 void
1551 NonlinearSystemBase::getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs)
1552 {
1553  const Node & node = _mesh.nodeRef(node_id);
1554  unsigned int s = number();
1555  if (node.has_dofs(s))
1556  {
1557  for (unsigned int v = 0; v < nVariables(); v++)
1558  for (unsigned int c = 0; c < node.n_comp(s, v); c++)
1559  dofs.push_back(node.dof_number(s, v, c));
1560  }
1561 }
1562 
1563 void
1565  GeometricSearchData & geom_search_data, std::map<dof_id_type, std::vector<dof_id_type>> & graph)
1566 {
1567  std::map<std::pair<unsigned int, unsigned int>, NearestNodeLocator *> & nearest_node_locators =
1568  geom_search_data._nearest_node_locators;
1569 
1570  const auto & node_to_elem_map = _mesh.nodeToElemMap();
1571  for (const auto & it : nearest_node_locators)
1572  {
1573  std::vector<dof_id_type> & slave_nodes = it.second->_slave_nodes;
1574 
1575  for (const auto & slave_node : slave_nodes)
1576  {
1577  std::set<dof_id_type> unique_slave_indices;
1578  std::set<dof_id_type> unique_master_indices;
1579 
1580  auto node_to_elem_pair = node_to_elem_map.find(slave_node);
1581  if (node_to_elem_pair != node_to_elem_map.end())
1582  {
1583  const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
1584 
1585  // Get the dof indices from each elem connected to the node
1586  for (const auto & cur_elem : elems)
1587  {
1588  std::vector<dof_id_type> dof_indices;
1589  dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
1590 
1591  for (const auto & dof : dof_indices)
1592  unique_slave_indices.insert(dof);
1593  }
1594  }
1595 
1596  std::vector<dof_id_type> master_nodes = it.second->_neighbor_nodes[slave_node];
1597 
1598  for (const auto & master_node : master_nodes)
1599  {
1600  auto master_node_to_elem_pair = node_to_elem_map.find(master_node);
1601  mooseAssert(master_node_to_elem_pair != node_to_elem_map.end(),
1602  "Missing entry in node to elem map");
1603  const std::vector<dof_id_type> & master_node_elems = master_node_to_elem_pair->second;
1604 
1605  // Get the dof indices from each elem connected to the node
1606  for (const auto & cur_elem : master_node_elems)
1607  {
1608  std::vector<dof_id_type> dof_indices;
1609  dofMap().dof_indices(_mesh.elemPtr(cur_elem), dof_indices);
1610 
1611  for (const auto & dof : dof_indices)
1612  unique_master_indices.insert(dof);
1613  }
1614  }
1615 
1616  for (const auto & slave_id : unique_slave_indices)
1617  for (const auto & master_id : unique_master_indices)
1618  {
1619  graph[slave_id].push_back(master_id);
1620  graph[master_id].push_back(slave_id);
1621  }
1622  }
1623  }
1624 
1625  // handle node-to-node constraints
1626  const auto & ncs = _constraints.getActiveNodalConstraints();
1627  for (const auto & nc : ncs)
1628  {
1629  std::vector<dof_id_type> master_dofs;
1630  std::vector<dof_id_type> & master_node_ids = nc->getMasterNodeId();
1631  for (const auto & node_id : master_node_ids)
1632  {
1633  Node * node = _mesh.queryNodePtr(node_id);
1634  if (node && node->processor_id() == this->processor_id())
1635  {
1636  getNodeDofs(node_id, master_dofs);
1637  }
1638  }
1639 
1640  _communicator.allgather(master_dofs);
1641 
1642  std::vector<dof_id_type> slave_dofs;
1643  std::vector<dof_id_type> & slave_node_ids = nc->getSlaveNodeId();
1644  for (const auto & node_id : slave_node_ids)
1645  {
1646  Node * node = _mesh.queryNodePtr(node_id);
1647  if (node && node->processor_id() == this->processor_id())
1648  {
1649  getNodeDofs(node_id, slave_dofs);
1650  }
1651  }
1652 
1653  _communicator.allgather(slave_dofs);
1654 
1655  for (const auto & master_id : master_dofs)
1656  for (const auto & slave_id : slave_dofs)
1657  {
1658  graph[master_id].push_back(slave_id);
1659  graph[slave_id].push_back(master_id);
1660  }
1661  }
1662 
1663  // Make every entry sorted and unique
1664  for (auto & it : graph)
1665  {
1666  std::vector<dof_id_type> & row = it.second;
1667  std::sort(row.begin(), row.end());
1668  std::vector<dof_id_type>::iterator uit = std::unique(row.begin(), row.end());
1669  row.resize(uit - row.begin());
1670  }
1671 }
1672 
1673 void
1675 {
1676  if (!hasMatrix(systemMatrixTag()))
1677  mooseError("Need a system matrix ");
1678 
1679  // At this point, have no idea how to make
1680  // this work with tag system
1681  auto & jacobian = getMatrix(systemMatrixTag());
1682 
1683  std::map<dof_id_type, std::vector<dof_id_type>> graph;
1684 
1685  findImplicitGeometricCouplingEntries(geom_search_data, graph);
1686 
1687  for (const auto & it : graph)
1688  {
1689  dof_id_type dof = it.first;
1690  const std::vector<dof_id_type> & row = it.second;
1691 
1692  for (const auto & coupled_dof : row)
1693  jacobian.add(dof, coupled_dof, 0);
1694  }
1695 }
1696 
1697 void
1699 {
1700  if (!hasMatrix(systemMatrixTag()))
1701  mooseError("A system matrix is required");
1702 
1703  auto & jacobian = getMatrix(systemMatrixTag());
1704 
1705 #if PETSC_VERSION_LESS_THAN(3, 3, 0)
1706 #else
1708  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
1709  MAT_NEW_NONZERO_ALLOCATION_ERR,
1710  PETSC_FALSE);
1712  MatSetOption(
1713  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
1714 #endif
1715 
1716  std::vector<numeric_index_type> zero_rows;
1717  std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators =
1718  NULL;
1719 
1720  if (!displaced)
1721  {
1722  GeometricSearchData & geom_search_data = _fe_problem.geomSearchData();
1723  penetration_locators = &geom_search_data._penetration_locators;
1724  }
1725  else
1726  {
1727  GeometricSearchData & displaced_geom_search_data =
1728  _fe_problem.getDisplacedProblem()->geomSearchData();
1729  penetration_locators = &displaced_geom_search_data._penetration_locators;
1730  }
1731 
1732  bool constraints_applied;
1734  constraints_applied = false;
1735  for (const auto & it : *penetration_locators)
1736  {
1738  {
1739  // Reset the constraint_applied flag before each new constraint, as they need to be
1740  // assembled separately
1741  constraints_applied = false;
1742  }
1743  PenetrationLocator & pen_loc = *(it.second);
1744 
1745  std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;
1746 
1747  BoundaryID slave_boundary = pen_loc._slave_boundary;
1748 
1749  zero_rows.clear();
1750  if (_constraints.hasActiveNodeFaceConstraints(slave_boundary, displaced))
1751  {
1752  const auto & constraints =
1753  _constraints.getActiveNodeFaceConstraints(slave_boundary, displaced);
1754 
1755  for (const auto & slave_node_num : slave_nodes)
1756  {
1757  Node & slave_node = _mesh.nodeRef(slave_node_num);
1758 
1759  if (slave_node.processor_id() == processor_id())
1760  {
1761  if (pen_loc._penetration_info[slave_node_num])
1762  {
1763  PenetrationInfo & info = *pen_loc._penetration_info[slave_node_num];
1764 
1765  const Elem * master_elem = info._elem;
1766  unsigned int master_side = info._side_num;
1767 
1768  // reinit variables at the node
1769  _fe_problem.reinitNodeFace(&slave_node, slave_boundary, 0);
1770 
1773 
1774  std::vector<Point> points;
1775  points.push_back(info._closest_point);
1776 
1777  // reinit variables on the master element's face at the contact point
1778  _fe_problem.setNeighborSubdomainID(master_elem, 0);
1779  _fe_problem.reinitNeighborPhys(master_elem, master_side, points, 0);
1780  for (const auto & nfc : constraints)
1781  {
1782  nfc->_jacobian = &jacobian;
1783 
1784  if (nfc->shouldApply())
1785  {
1786  constraints_applied = true;
1787 
1788  nfc->subProblem().prepareShapes(nfc->variable().number(), 0);
1789  nfc->subProblem().prepareNeighborShapes(nfc->variable().number(), 0);
1790 
1791  nfc->computeJacobian();
1792 
1793  if (nfc->overwriteSlaveJacobian())
1794  {
1795  // Add this variable's dof's row to be zeroed
1796  zero_rows.push_back(nfc->variable().nodalDofIndex());
1797  }
1798 
1799  std::vector<dof_id_type> slave_dofs(1, nfc->variable().nodalDofIndex());
1800 
1801  // Cache the jacobian block for the slave side
1803  slave_dofs,
1804  nfc->_connected_dof_indices,
1805  nfc->variable().scalingFactor());
1806 
1807  // Cache the jacobian block for the master side
1808  if (nfc->addCouplingEntriesToJacobian())
1810  nfc->_Kne,
1811  nfc->masterVariable().dofIndicesNeighbor(),
1812  nfc->_connected_dof_indices,
1813  nfc->variable().scalingFactor());
1814 
1816  if (nfc->addCouplingEntriesToJacobian())
1818 
1819  // Do the off-diagonals next
1820  const std::vector<MooseVariableFEBase *> coupled_vars = nfc->getCoupledMooseVars();
1821  for (const auto & jvar : coupled_vars)
1822  {
1823  // Only compute jacobians for nonlinear variables
1824  if (jvar->kind() != Moose::VAR_NONLINEAR)
1825  continue;
1826 
1827  // Only compute Jacobian entries if this coupling is being used by the
1828  // preconditioner
1829  if (nfc->variable().number() == jvar->number() ||
1830  !_fe_problem.areCoupled(nfc->variable().number(), jvar->number()))
1831  continue;
1832 
1833  // Need to zero out the matrices first
1835 
1836  nfc->subProblem().prepareShapes(nfc->variable().number(), 0);
1837  nfc->subProblem().prepareNeighborShapes(jvar->number(), 0);
1838 
1839  nfc->computeOffDiagJacobian(jvar->number());
1840 
1841  // Cache the jacobian block for the slave side
1843  slave_dofs,
1844  nfc->_connected_dof_indices,
1845  nfc->variable().scalingFactor());
1846 
1847  // Cache the jacobian block for the master side
1848  if (nfc->addCouplingEntriesToJacobian())
1850  nfc->variable().dofIndicesNeighbor(),
1851  nfc->_connected_dof_indices,
1852  nfc->variable().scalingFactor());
1853 
1855  if (nfc->addCouplingEntriesToJacobian())
1857  }
1858  }
1859  }
1860  }
1861  }
1862  }
1863  }
1865  {
1866  // See if constraints were applied anywhere
1867  _communicator.max(constraints_applied);
1868 
1869  if (constraints_applied)
1870  {
1871 #ifdef LIBMESH_HAVE_PETSC
1872 // Necessary for speed
1873 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
1874  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS);
1875 #elif PETSC_VERSION_LESS_THAN(3, 1, 0)
1876  // In Petsc 3.0.0, MatSetOption has three args...the third arg
1877  // determines whether the option is set (true) or unset (false)
1878  MatSetOption(
1879  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
1880 #else
1881  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
1882  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
1883  PETSC_TRUE);
1884 #endif
1885 #endif
1886 
1887  jacobian.close();
1888  jacobian.zero_rows(zero_rows, 0.0);
1889  jacobian.close();
1891  jacobian.close();
1892  }
1893  }
1894  }
1896  {
1897  // See if constraints were applied anywhere
1898  _communicator.max(constraints_applied);
1899 
1900  if (constraints_applied)
1901  {
1902 #ifdef LIBMESH_HAVE_PETSC
1903 // Necessary for speed
1904 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
1905  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS);
1906 #elif PETSC_VERSION_LESS_THAN(3, 1, 0)
1907  // In Petsc 3.0.0, MatSetOption has three args...the third arg
1908  // determines whether the option is set (true) or unset (false)
1909  MatSetOption(
1910  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
1911 #else
1912  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
1913  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
1914  PETSC_TRUE);
1915 #endif
1916 #endif
1917 
1918  jacobian.close();
1919  jacobian.zero_rows(zero_rows, 0.0);
1920  jacobian.close();
1922  jacobian.close();
1923  }
1924  }
1925 
1926  mortarJacobianConstraints(displaced);
1927 
1928  THREAD_ID tid = 0;
1929  // go over element-element constraint interface
1930  std::map<unsigned int, std::shared_ptr<ElementPairLocator>> * element_pair_locators = nullptr;
1931 
1932  if (!displaced)
1933  {
1934  GeometricSearchData & geom_search_data = _fe_problem.geomSearchData();
1935  element_pair_locators = &geom_search_data._element_pair_locators;
1936  }
1937  else
1938  {
1939  GeometricSearchData & displaced_geom_search_data =
1940  _fe_problem.getDisplacedProblem()->geomSearchData();
1941  element_pair_locators = &displaced_geom_search_data._element_pair_locators;
1942  }
1943 
1944  for (const auto & it : *element_pair_locators)
1945  {
1946  ElementPairLocator & elem_pair_loc = *(it.second);
1947 
1948  if (_constraints.hasActiveElemElemConstraints(it.first, displaced))
1949  {
1950  // ElemElemConstraint objects
1951  const auto & _element_constraints =
1952  _constraints.getActiveElemElemConstraints(it.first, displaced);
1953 
1954  // go over pair elements
1955  const std::list<std::pair<const Elem *, const Elem *>> & elem_pairs =
1956  elem_pair_loc.getElemPairs();
1957  for (const auto & pr : elem_pairs)
1958  {
1959  const Elem * elem1 = pr.first;
1960  const Elem * elem2 = pr.second;
1961 
1962  if (elem1->processor_id() != processor_id())
1963  continue;
1964 
1965  const ElementPairInfo & info = elem_pair_loc.getElemPairInfo(pr);
1966 
1967  // for each element process constraints on the
1968  for (const auto & ec : _element_constraints)
1969  {
1970  _fe_problem.setCurrentSubdomainID(elem1, tid);
1972  _fe_problem.setNeighborSubdomainID(elem2, tid);
1974 
1975  ec->subProblem().prepareShapes(ec->variable().number(), tid);
1976  ec->subProblem().prepareNeighborShapes(ec->variable().number(), tid);
1977 
1978  ec->reinit(info);
1979  ec->computeJacobian();
1982  }
1984  }
1985  }
1986  }
1987 
1988  // go over NodeELemConstraints
1989  std::set<dof_id_type> unique_slave_node_ids;
1990  constraints_applied = false;
1991  for (const auto & slave_id : _mesh.meshSubdomains())
1992  {
1993  for (const auto & master_id : _mesh.meshSubdomains())
1994  {
1995  if (_constraints.hasActiveNodeElemConstraints(slave_id, master_id, displaced))
1996  {
1997  const auto & constraints =
1998  _constraints.getActiveNodeElemConstraints(slave_id, master_id, displaced);
1999 
2000  // get unique set of ids of all nodes on current block
2001  unique_slave_node_ids.clear();
2002  const MeshBase & meshhelper = _mesh.getMesh();
2003  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(slave_id),
2004  meshhelper.active_subdomain_elements_end(slave_id)))
2005  {
2006  for (auto & n : elem->node_ref_range())
2007  unique_slave_node_ids.insert(n.id());
2008  }
2009 
2010  for (auto slave_node_id : unique_slave_node_ids)
2011  {
2012  const Node & slave_node = _mesh.nodeRef(slave_node_id);
2013  // check if slave node is on current processor
2014  if (slave_node.processor_id() == processor_id())
2015  {
2016  // This reinits the variables that exist on the slave node
2017  _fe_problem.reinitNodeFace(&slave_node, slave_id, 0);
2018 
2019  // This will set aside residual and jacobian space for the variables that have dofs
2020  // on the slave node
2023 
2024  for (const auto & nec : constraints)
2025  {
2026  if (nec->shouldApply())
2027  {
2028  constraints_applied = true;
2029 
2030  nec->_jacobian = &jacobian;
2031  nec->subProblem().prepareShapes(nec->variable().number(), 0);
2032  nec->subProblem().prepareNeighborShapes(nec->variable().number(), 0);
2033 
2034  nec->computeJacobian();
2035 
2036  if (nec->overwriteSlaveJacobian())
2037  {
2038  // Add this variable's dof's row to be zeroed
2039  zero_rows.push_back(nec->variable().nodalDofIndex());
2040  }
2041 
2042  std::vector<dof_id_type> slave_dofs(1, nec->variable().nodalDofIndex());
2043 
2044  // Cache the jacobian block for the slave side
2046  slave_dofs,
2047  nec->_connected_dof_indices,
2048  nec->variable().scalingFactor());
2049 
2050  // Cache the jacobian block for the master side
2052  nec->_Kne,
2053  nec->masterVariable().dofIndicesNeighbor(),
2054  nec->_connected_dof_indices,
2055  nec->variable().scalingFactor());
2056 
2059 
2060  // Do the off-diagonals next
2061  const std::vector<MooseVariableFEBase *> coupled_vars = nec->getCoupledMooseVars();
2062  for (const auto & jvar : coupled_vars)
2063  {
2064  // Only compute jacobians for nonlinear variables
2065  if (jvar->kind() != Moose::VAR_NONLINEAR)
2066  continue;
2067 
2068  // Only compute Jacobian entries if this coupling is being used by the
2069  // preconditioner
2070  if (nec->variable().number() == jvar->number() ||
2071  !_fe_problem.areCoupled(nec->variable().number(), jvar->number()))
2072  continue;
2073 
2074  // Need to zero out the matrices first
2076 
2077  nec->subProblem().prepareShapes(nec->variable().number(), 0);
2078  nec->subProblem().prepareNeighborShapes(jvar->number(), 0);
2079 
2080  nec->computeOffDiagJacobian(jvar->number());
2081 
2082  // Cache the jacobian block for the slave side
2084  slave_dofs,
2085  nec->_connected_dof_indices,
2086  nec->variable().scalingFactor());
2087 
2088  // Cache the jacobian block for the master side
2090  nec->variable().dofIndicesNeighbor(),
2091  nec->_connected_dof_indices,
2092  nec->variable().scalingFactor());
2093 
2096  }
2097  }
2098  }
2099  }
2100  }
2101  }
2102  }
2103  }
2104  // See if constraints were applied anywhere
2105  _communicator.max(constraints_applied);
2106 
2107  if (constraints_applied)
2108  {
2109 #ifdef LIBMESH_HAVE_PETSC
2110 // Necessary for speed
2111 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
2112  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS);
2113 #elif PETSC_VERSION_LESS_THAN(3, 1, 0)
2114  // In Petsc 3.0.0, MatSetOption has three args...the third arg
2115  // determines whether the option is set (true) or unset (false)
2116  MatSetOption(
2117  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
2118 #else
2119  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2120  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2121  PETSC_TRUE);
2122 #endif
2123 #endif
2124 
2125  jacobian.close();
2126  jacobian.zero_rows(zero_rows, 0.0);
2127  jacobian.close();
2129  jacobian.close();
2130  }
2131 }
2132 
2133 void
2135 {
2136  MooseObjectWarehouse<ScalarKernel> * scalar_kernel_warehouse;
2137 
2138  if (!tags.size() || tags.size() == _fe_problem.numMatrixTags())
2139  scalar_kernel_warehouse = &_scalar_kernels;
2140  else if (tags.size() == 1)
2141  scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2142  else
2143  scalar_kernel_warehouse = &(_scalar_kernels.getMatrixTagsObjectWarehouse(tags, 0));
2144 
2145  // Compute the diagonal block for scalar variables
2146  if (scalar_kernel_warehouse->hasActiveObjects())
2147  {
2148  const auto & scalars = scalar_kernel_warehouse->getActiveObjects();
2149 
2150  _fe_problem.reinitScalars(/*tid=*/0);
2151 
2152  bool have_scalar_contributions = false;
2153  for (const auto & kernel : scalars)
2154  {
2155  kernel->reinit();
2156  const std::vector<dof_id_type> & dof_indices = kernel->variable().dofIndices();
2157  const DofMap & dof_map = kernel->variable().dofMap();
2158  const dof_id_type first_dof = dof_map.first_dof();
2159  const dof_id_type end_dof = dof_map.end_dof();
2160  for (dof_id_type dof : dof_indices)
2161  {
2162  if (dof >= first_dof && dof < end_dof)
2163  {
2164  kernel->computeJacobian();
2165  _fe_problem.addJacobianOffDiagScalar(kernel->variable().number());
2166  have_scalar_contributions = true;
2167  break;
2168  }
2169  }
2170  }
2171 
2172  if (have_scalar_contributions)
2174  }
2175 }
2176 
2177 void
2179 {
2180  // Make matrix ready to use
2182 
2183  for (auto tag : tags)
2184  {
2185  if (!hasMatrix(tag))
2186  continue;
2187 
2188  auto & jacobian = getMatrix(tag);
2189 #ifdef LIBMESH_HAVE_PETSC
2190 // Necessary for speed
2191 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
2192  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS);
2193 #elif PETSC_VERSION_LESS_THAN(3, 1, 0)
2194  // In Petsc 3.0.0, MatSetOption has three args...the third arg
2195  // determines whether the option is set (true) or unset (false)
2196  MatSetOption(
2197  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
2198 #else
2199  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2200  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2201  PETSC_TRUE);
2202 #endif
2203 #if PETSC_VERSION_LESS_THAN(3, 3, 0)
2204 #else
2206  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2207  MAT_NEW_NONZERO_ALLOCATION_ERR,
2208  PETSC_FALSE);
2209 #endif
2210 
2211 #endif
2212  }
2213  // jacobianSetup /////
2214  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2215  {
2216  _kernels.jacobianSetup(tid);
2220  if (_doing_dg)
2226  }
2231 
2232  // reinit scalar variables
2233  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2235 
2236  PARALLEL_TRY
2237  {
2238  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
2239  switch (_fe_problem.coupling())
2240  {
2241  case Moose::COUPLING_DIAG:
2242  {
2244  Threads::parallel_reduce(elem_range, cj);
2245 
2246  unsigned int n_threads = libMesh::n_threads();
2247  for (unsigned int i = 0; i < n_threads;
2248  i++) // Add any Jacobian contributions still hanging around
2250 
2251  // Block restricted Nodal Kernels
2253  {
2255  ConstNodeRange & range = *_mesh.getLocalNodeRange();
2256  Threads::parallel_reduce(range, cnkjt);
2257 
2258  unsigned int n_threads = libMesh::n_threads();
2259  for (unsigned int i = 0; i < n_threads;
2260  i++) // Add any cached jacobians that might be hanging around
2262  }
2263 
2264  // Boundary restricted Nodal Kernels
2266  {
2268  ConstBndNodeRange & bnd_range = *_mesh.getBoundaryNodeRange();
2269 
2270  Threads::parallel_reduce(bnd_range, cnkjt);
2271  unsigned int n_threads = libMesh::n_threads();
2272  for (unsigned int i = 0; i < n_threads;
2273  i++) // Add any cached jacobians that might be hanging around
2275  }
2276  }
2277  break;
2278 
2279  default:
2281  {
2283  Threads::parallel_reduce(elem_range, cj);
2284  unsigned int n_threads = libMesh::n_threads();
2285 
2286  for (unsigned int i = 0; i < n_threads; i++)
2288 
2289  // Block restricted Nodal Kernels
2291  {
2293  ConstNodeRange & range = *_mesh.getLocalNodeRange();
2294  Threads::parallel_reduce(range, cnkjt);
2295 
2296  unsigned int n_threads = libMesh::n_threads();
2297  for (unsigned int i = 0; i < n_threads;
2298  i++) // Add any cached jacobians that might be hanging around
2300  }
2301 
2302  // Boundary restricted Nodal Kernels
2304  {
2306  ConstBndNodeRange & bnd_range = *_mesh.getBoundaryNodeRange();
2307 
2308  Threads::parallel_reduce(bnd_range, cnkjt);
2309 
2310  unsigned int n_threads = libMesh::n_threads();
2311  for (unsigned int i = 0; i < n_threads;
2312  i++) // Add any cached jacobians that might be hanging around
2314  }
2315  }
2316  break;
2317  }
2318 
2319  computeDiracContributions(tags, true);
2320 
2322 
2323  static bool first = true;
2324 
2325  // This adds zeroes into geometric coupling entries to ensure they stay in the matrix
2327  {
2328  first = false;
2330 
2333  }
2334  }
2335  PARALLEL_CATCH;
2336 
2337  closeTaggedMatrices(tags);
2338 
2339  // Have no idea how to have constraints work
2340  // with the tag system
2341  PARALLEL_TRY
2342  {
2343  // Add in Jacobian contributions from Constraints
2345  {
2346  // Nodal Constraints
2348 
2349  // Undisplaced Constraints
2350  constraintJacobians(false);
2351 
2352  // Displaced Constraints
2354  constraintJacobians(true);
2355  }
2356  }
2357  PARALLEL_CATCH;
2359  closeTaggedMatrices(tags);
2360 
2361  // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
2362  // on boundary nodes
2363  if (_has_diag_save_in)
2365 
2366  PARALLEL_TRY
2367  {
2368  MooseObjectWarehouse<NodalBCBase> * nbc_warehouse;
2369  // Select nodal kernels
2370  if (tags.size() == _fe_problem.numMatrixTags() || !tags.size())
2371  nbc_warehouse = &_nodal_bcs;
2372  else if (tags.size() == 1)
2373  nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0));
2374  else
2375  nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0));
2376 
2377  // Cache the information about which BCs are coupled to which
2378  // variables, so we don't have to figure it out for each node.
2379  std::map<std::string, std::set<unsigned int>> bc_involved_vars;
2380  const std::set<BoundaryID> & all_boundary_ids = _mesh.getBoundaryIDs();
2381  for (const auto & bid : all_boundary_ids)
2382  {
2383  // Get reference to all the NodalBCs for this ID. This is only
2384  // safe if there are NodalBCBases there to be gotten...
2385  if (nbc_warehouse->hasActiveBoundaryObjects(bid))
2386  {
2387  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid);
2388  for (const auto & bc : bcs)
2389  {
2390  const std::vector<MooseVariableFEBase *> & coupled_moose_vars = bc->getCoupledMooseVars();
2391 
2392  // Create the set of "involved" MOOSE nonlinear vars, which includes all coupled vars
2393  // and the BC's own variable
2394  std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
2395  for (const auto & coupled_var : coupled_moose_vars)
2396  if (coupled_var->kind() == Moose::VAR_NONLINEAR)
2397  var_set.insert(coupled_var->number());
2398 
2399  var_set.insert(bc->variable().number());
2400  }
2401  }
2402  }
2403 
2404  // Get variable coupling list. We do all the NodalBCBase stuff on
2405  // thread 0... The couplingEntries() data structure determines
2406  // which variables are "coupled" as far as the preconditioner is
2407  // concerned, not what variables a boundary condition specifically
2408  // depends on.
2409  std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> & coupling_entries =
2410  _fe_problem.couplingEntries(/*_tid=*/0);
2411 
2412  // Compute Jacobians for NodalBCBases
2413  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
2414  for (const auto & bnode : bnd_nodes)
2415  {
2416  BoundaryID boundary_id = bnode->_bnd_id;
2417  Node * node = bnode->_node;
2418 
2419  if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) &&
2420  node->processor_id() == processor_id())
2421  {
2422  _fe_problem.reinitNodeFace(node, boundary_id, 0);
2423 
2424  const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id);
2425  for (const auto & bc : bcs)
2426  {
2427  // Get the set of involved MOOSE vars for this BC
2428  std::set<unsigned int> & var_set = bc_involved_vars[bc->name()];
2429 
2430  // Loop over all the variables whose Jacobian blocks are
2431  // actually being computed, call computeOffDiagJacobian()
2432  // for each one which is actually coupled (otherwise the
2433  // value is zero.)
2434  for (const auto & it : coupling_entries)
2435  {
2436  unsigned int ivar = it.first->number(), jvar = it.second->number();
2437 
2438  // We are only going to call computeOffDiagJacobian() if:
2439  // 1.) the BC's variable is ivar
2440  // 2.) jvar is "involved" with the BC (including jvar==ivar), and
2441  // 3.) the BC should apply.
2442  if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply())
2443  bc->computeOffDiagJacobian(jvar);
2444  }
2445  }
2446  }
2447  } // end loop over boundary nodes
2448 
2449  // Set the cached NodalBCBase values in the Jacobian matrix
2451  }
2452  PARALLEL_CATCH;
2453 
2454  closeTaggedMatrices(tags);
2455 
2456  // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs
2457  // on boundary nodes
2460 
2461  if (hasDiagSaveIn())
2463 }
2464 
2465 void
2466 NonlinearSystemBase::setVariableGlobalDoFs(const std::string & var_name)
2467 {
2468  AllLocalDofIndicesThread aldit(_sys.system(), {var_name});
2469  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
2470  Threads::parallel_reduce(elem_range, aldit);
2471  _communicator.set_union(aldit._all_dof_indices);
2472  _var_all_dof_indices.assign(aldit._all_dof_indices.begin(), aldit._all_dof_indices.end());
2473 }
2474 
2475 void
2476 NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian)
2477 {
2478  _nl_matrix_tags.clear();
2479 
2480  auto & tags = _fe_problem.getMatrixTags();
2481 
2482  for (auto & tag : tags)
2483  _nl_matrix_tags.insert(tag.second);
2484 
2485  computeJacobian(jacobian, _nl_matrix_tags);
2486 }
2487 
2488 void
2489 NonlinearSystemBase::computeJacobian(SparseMatrix<Number> & jacobian, const std::set<TagID> & tags)
2490 {
2492 
2493  computeJacobianTags(tags);
2494 
2496 }
2497 
2498 void
2499 NonlinearSystemBase::computeJacobianTags(const std::set<TagID> & tags)
2500 {
2501  TIME_SECTION(_compute_jacobian_tags_timer);
2502 
2503  FloatingPointExceptionGuard fpe_guard(_app);
2504 
2505  try
2506  {
2508  }
2509  catch (MooseException & e)
2510  {
2511  // The buck stops here, we have already handled the exception by
2512  // calling stopSolve(), it is now up to PETSc to return a
2513  // "diverged" reason during the next solve.
2514  }
2515 }
2516 
2517 void
2518 NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
2519 {
2520  _nl_matrix_tags.clear();
2521 
2522  auto & tags = _fe_problem.getMatrixTags();
2523  for (auto & tag : tags)
2524  _nl_matrix_tags.insert(tag.second);
2525 
2527 }
2528 
2529 void
2530 NonlinearSystemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
2531  const std::set<TagID> & tags)
2532 {
2533  TIME_SECTION(_compute_jacobian_blocks_timer);
2534 
2535  FloatingPointExceptionGuard fpe_guard(_app);
2536 
2537  for (unsigned int i = 0; i < blocks.size(); i++)
2538  {
2539  SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
2540 
2541 #ifdef LIBMESH_HAVE_PETSC
2542 // Necessary for speed
2543 #if PETSC_VERSION_LESS_THAN(3, 0, 0)
2544  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS);
2545 #elif PETSC_VERSION_LESS_THAN(3, 1, 0)
2546  // In Petsc 3.0.0, MatSetOption has three args...the third arg
2547  // determines whether the option is set (true) or unset (false)
2548  MatSetOption(
2549  static_cast<PetscMatrix<Number> &>(jacobian).mat(), MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
2550 #else
2551  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2552  MAT_KEEP_NONZERO_PATTERN, // This is changed in 3.1
2553  PETSC_TRUE);
2554 #endif
2555 #if PETSC_VERSION_LESS_THAN(3, 3, 0)
2556 #else
2558  MatSetOption(static_cast<PetscMatrix<Number> &>(jacobian).mat(),
2559  MAT_NEW_NONZERO_ALLOCATION_ERR,
2560  PETSC_FALSE);
2561 #endif
2562 
2563 #endif
2564 
2565  jacobian.zero();
2566  }
2567 
2568  for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++)
2570 
2571  PARALLEL_TRY
2572  {
2573  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
2574  ComputeJacobianBlocksThread cjb(_fe_problem, blocks, tags);
2575  Threads::parallel_reduce(elem_range, cjb);
2576  }
2577  PARALLEL_CATCH;
2578 
2579  for (unsigned int i = 0; i < blocks.size(); i++)
2580  blocks[i]->_jacobian.close();
2581 
2582  for (unsigned int i = 0; i < blocks.size(); i++)
2583  {
2584  libMesh::System & precond_system = blocks[i]->_precond_system;
2585  SparseMatrix<Number> & jacobian = blocks[i]->_jacobian;
2586 
2587  unsigned int ivar = blocks[i]->_ivar;
2588  unsigned int jvar = blocks[i]->_jvar;
2589 
2590  // Dirichlet BCs
2591  std::vector<numeric_index_type> zero_rows;
2592  PARALLEL_TRY
2593  {
2594  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
2595  for (const auto & bnode : bnd_nodes)
2596  {
2597  BoundaryID boundary_id = bnode->_bnd_id;
2598  Node * node = bnode->_node;
2599 
2600  if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
2601  {
2602  const auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
2603 
2604  if (node->processor_id() == processor_id())
2605  {
2606  _fe_problem.reinitNodeFace(node, boundary_id, 0);
2607 
2608  for (const auto & bc : bcs)
2609  if (bc->variable().number() == ivar && bc->shouldApply())
2610  {
2611  // The first zero is for the variable number... there is only one variable in each
2612  // mini-system
2613  // The second zero only works with Lagrange elements!
2614  zero_rows.push_back(node->dof_number(precond_system.number(), 0, 0));
2615  }
2616  }
2617  }
2618  }
2619  }
2620  PARALLEL_CATCH;
2621 
2622  jacobian.close();
2623 
2624  // This zeroes the rows corresponding to Dirichlet BCs and puts a 1.0 on the diagonal
2625  if (ivar == jvar)
2626  jacobian.zero_rows(zero_rows, 1.0);
2627  else
2628  jacobian.zero_rows(zero_rows, 0.0);
2629 
2630  jacobian.close();
2631  }
2632 }
2633 
2634 void
2636 {
2643  _kernels.updateActive(tid);
2646  if (tid == 0)
2647  {
2654  }
2655 }
2656 
2657 Real
2658 NonlinearSystemBase::computeDamping(const NumericVector<Number> & solution,
2659  const NumericVector<Number> & update)
2660 {
2661  // Default to no damping
2662  Real damping = 1.0;
2663  bool has_active_dampers = false;
2664 
2666  {
2667  TIME_SECTION(_compute_dampers_timer);
2668 
2669  has_active_dampers = true;
2672  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cid);
2673  damping = std::min(cid.damping(), damping);
2674  }
2675 
2677  {
2678  TIME_SECTION(_compute_dampers_timer);
2679 
2680  has_active_dampers = true;
2683  Threads::parallel_reduce(*_mesh.getLocalNodeRange(), cndt);
2684  damping = std::min(cndt.damping(), damping);
2685  }
2686 
2688  {
2689  TIME_SECTION(_compute_dampers_timer);
2690 
2691  has_active_dampers = true;
2692  const auto & gdampers = _general_dampers.getActiveObjects();
2693  for (const auto & damper : gdampers)
2694  {
2695  Real gd_damping = damper->computeDamping(solution, update);
2696  try
2697  {
2698  damper->checkMinDamping(gd_damping);
2699  }
2700  catch (MooseException & e)
2701  {
2703  }
2704  damping = std::min(gd_damping, damping);
2705  }
2706  }
2707 
2708  _communicator.min(damping);
2709 
2710  if (has_active_dampers && damping < 1.0)
2711  _console << " Damping factor: " << damping << "\n";
2712 
2713  return damping;
2714 }
2715 
2716 void
2717 NonlinearSystemBase::computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian)
2718 {
2720 
2721  std::set<const Elem *> dirac_elements;
2722 
2724  {
2725  TIME_SECTION(_compute_dirac_timer);
2726 
2727  // TODO: Need a threading fix... but it's complicated!
2728  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2729  {
2730  const auto & dkernels = _dirac_kernels.getActiveObjects(tid);
2731  for (const auto & dkernel : dkernels)
2732  {
2733  dkernel->clearPoints();
2734  dkernel->addPoints();
2735  }
2736  }
2737 
2738  ComputeDiracThread cd(_fe_problem, tags, is_jacobian);
2739 
2740  _fe_problem.getDiracElements(dirac_elements);
2741 
2742  DistElemRange range(dirac_elements.begin(), dirac_elements.end(), 1);
2743  // TODO: Make Dirac work thread!
2744  // Threads::parallel_reduce(range, cd);
2745 
2746  cd(range);
2747  }
2748 }
2749 
2750 NumericVector<Number> &
2752 {
2753  _need_residual_copy = true;
2754  return _residual_copy;
2755 }
2756 
2757 NumericVector<Number> &
2759 {
2760  _need_residual_ghosted = true;
2761  if (!_residual_ghosted)
2762  _residual_ghosted = &addVector("residual_ghosted", false, GHOSTED);
2763  return *_residual_ghosted;
2764 }
2765 
2766 void
2767 NonlinearSystemBase::augmentSparsity(SparsityPattern::Graph & sparsity,
2768  std::vector<dof_id_type> & n_nz,
2769  std::vector<dof_id_type> & n_oz)
2770 {
2772  {
2774 
2775  std::map<dof_id_type, std::vector<dof_id_type>> graph;
2776 
2778 
2781  graph);
2782 
2783  const dof_id_type first_dof_on_proc = dofMap().first_dof(processor_id());
2784  const dof_id_type end_dof_on_proc = dofMap().end_dof(processor_id());
2785 
2786  // The total number of dofs on and off processor
2787  const dof_id_type n_dofs_on_proc = dofMap().n_local_dofs();
2788  const dof_id_type n_dofs_not_on_proc = dofMap().n_dofs() - dofMap().n_local_dofs();
2789 
2790  for (const auto & git : graph)
2791  {
2792  dof_id_type dof = git.first;
2793  dof_id_type local_dof = dof - first_dof_on_proc;
2794 
2795  if (dof < first_dof_on_proc || dof >= end_dof_on_proc)
2796  continue;
2797 
2798  const std::vector<dof_id_type> & row = git.second;
2799 
2800  SparsityPattern::Row & sparsity_row = sparsity[local_dof];
2801 
2802  unsigned int original_row_length = sparsity_row.size();
2803 
2804  sparsity_row.insert(sparsity_row.end(), row.begin(), row.end());
2805 
2806  SparsityPattern::sort_row(
2807  sparsity_row.begin(), sparsity_row.begin() + original_row_length, sparsity_row.end());
2808 
2809  // Fix up nonzero arrays
2810  for (const auto & coupled_dof : row)
2811  {
2812  if (coupled_dof < first_dof_on_proc || coupled_dof >= end_dof_on_proc)
2813  {
2814  if (n_oz[local_dof] < n_dofs_not_on_proc)
2815  n_oz[local_dof]++;
2816  }
2817  else
2818  {
2819  if (n_nz[local_dof] < n_dofs_on_proc)
2820  n_nz[local_dof]++;
2821  }
2822  }
2823  }
2824  }
2825 }
2826 
2827 void
2829 {
2831  {
2832  if (!_serialized_solution.initialized() || _serialized_solution.size() != _sys.n_dofs())
2833  {
2834  _serialized_solution.clear();
2835  _serialized_solution.init(_sys.n_dofs(), false, SERIAL);
2836  }
2837 
2839  }
2840 }
2841 
2842 void
2843 NonlinearSystemBase::setSolution(const NumericVector<Number> & soln)
2844 {
2845  _current_solution = &soln;
2846 
2849 }
2850 
2851 void
2852 NonlinearSystemBase::setSolutionUDot(const NumericVector<Number> & u_dot)
2853 {
2854  *_u_dot = u_dot;
2855 }
2856 
2857 void
2858 NonlinearSystemBase::setSolutionUDotDot(const NumericVector<Number> & u_dotdot)
2859 {
2860  *_u_dotdot = u_dotdot;
2861 }
2862 
2863 void
2864 NonlinearSystemBase::setSolutionUDotOld(const NumericVector<Number> & u_dot_old)
2865 {
2866  *_u_dot_old = u_dot_old;
2867 }
2868 
2869 void
2870 NonlinearSystemBase::setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old)
2871 {
2872  *_u_dotdot_old = u_dotdot_old;
2873 }
2874 
2875 NumericVector<Number> &
2877 {
2878  if (!_serialized_solution.initialized())
2879  _serialized_solution.init(_sys.n_dofs(), false, SERIAL);
2880 
2882  return _serialized_solution;
2883 }
2884 
2885 void
2886 NonlinearSystemBase::setPreconditioner(std::shared_ptr<MoosePreconditioner> pc)
2887 {
2888  if (_preconditioner.get() != nullptr)
2889  mooseError("More than one active Preconditioner detected");
2890 
2891  _preconditioner = pc;
2892 }
2893 
2894 MoosePreconditioner const *
2896 {
2897  return _preconditioner.get();
2898 }
2899 
2900 void
2902 {
2903  _increment_vec = &_sys.add_vector("u_increment", true, GHOSTED);
2904 }
2905 
2906 void
2908  const std::set<MooseVariable *> & damped_vars)
2909 {
2910  for (const auto & var : damped_vars)
2911  var->computeIncrementAtQps(*_increment_vec);
2912 }
2913 
2914 void
2916  const std::set<MooseVariable *> & damped_vars)
2917 {
2918  for (const auto & var : damped_vars)
2919  var->computeIncrementAtNode(*_increment_vec);
2920 }
2921 
2922 void
2923 NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const
2924 {
2925  // Check kernel coverage of subdomains (blocks) in your mesh
2926  std::set<SubdomainID> input_subdomains;
2927  std::set<std::string> kernel_variables;
2928 
2929  bool global_kernels_exist = _kernels.hasActiveBlockObjects(Moose::ANY_BLOCK_ID);
2930  global_kernels_exist |= _scalar_kernels.hasActiveObjects();
2931  global_kernels_exist |= _nodal_kernels.hasActiveObjects();
2932 
2933  _kernels.subdomainsCovered(input_subdomains, kernel_variables);
2934  _nodal_kernels.subdomainsCovered(input_subdomains, kernel_variables);
2935  _scalar_kernels.subdomainsCovered(input_subdomains, kernel_variables);
2936  _constraints.subdomainsCovered(input_subdomains, kernel_variables);
2937 
2938  if (!global_kernels_exist)
2939  {
2940  std::set<SubdomainID> difference;
2941  std::set_difference(mesh_subdomains.begin(),
2942  mesh_subdomains.end(),
2943  input_subdomains.begin(),
2944  input_subdomains.end(),
2945  std::inserter(difference, difference.end()));
2946 
2947  if (!difference.empty())
2948  {
2949  std::stringstream missing_block_ids;
2950  std::copy(difference.begin(),
2951  difference.end(),
2952  std::ostream_iterator<unsigned int>(missing_block_ids, " "));
2953  mooseError("Each subdomain must contain at least one Kernel.\nThe following block(s) lack an "
2954  "active kernel: " +
2955  missing_block_ids.str());
2956  }
2957  }
2958 
2959  std::set<VariableName> variables(getVariableNames().begin(), getVariableNames().end());
2960 
2961  std::set<VariableName> difference;
2962  std::set_difference(variables.begin(),
2963  variables.end(),
2964  kernel_variables.begin(),
2965  kernel_variables.end(),
2966  std::inserter(difference, difference.end()));
2967 
2968  if (!difference.empty())
2969  {
2970  std::stringstream missing_kernel_vars;
2971  std::copy(difference.begin(),
2972  difference.end(),
2973  std::ostream_iterator<std::string>(missing_kernel_vars, " "));
2974  mooseError("Each variable must be referenced by at least one active Kernel.\nThe following "
2975  "variable(s) lack an active kernel: " +
2976  missing_kernel_vars.str());
2977  }
2978 }
2979 
2980 bool
2982 {
2983  auto & time_kernels = _kernels.getVectorTagObjectWarehouse(timeVectorTag(), 0);
2984 
2985  return time_kernels.hasActiveObjects();
2986 }
2987 
2988 void
2990 {
2991  if (pcs == "left")
2993  else if (pcs == "right")
2995  else if (pcs == "symmetric")
2997  else if (pcs == "default")
2999  else
3000  mooseError("Unknown PC side specified.");
3001 }
3002 
3003 void
3005 {
3006  if (kspnorm == "none")
3008  else if (kspnorm == "preconditioned")
3010  else if (kspnorm == "unpreconditioned")
3012  else if (kspnorm == "natural")
3014  else if (kspnorm == "default")
3016  else
3017  mooseError("Unknown ksp norm type specified.");
3018 }
3019 
3020 bool
3022 {
3023  return _integrated_bcs.hasActiveBoundaryObjects(bnd_id, tid);
3024 }
3025 
3027  THREAD_ID /*tid*/) const
3028 {
3029  return _doing_dg;
3030 }
3031 
3032 bool
3034 {
3035  return _doing_dg;
3036 }
3037 
3038 void
3039 NonlinearSystemBase::setPreviousNewtonSolution(const NumericVector<Number> & soln)
3040 {
3042  *_solution_previous_nl = soln;
3043 }
3044 
3045 void
3047 {
3048  // go over mortar constraints
3049  const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
3050 
3051  std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor<RESIDUAL>>::iterator
3052  it,
3053  end_it;
3054 
3055  for (const auto & mortar_interface : mortar_interfaces)
3056  {
3057  if (!displaced)
3058  {
3059  it = _undisplaced_mortar_residual_functors.find(mortar_interface.first);
3061  }
3062  else
3063  {
3064  it = _displaced_mortar_residual_functors.find(mortar_interface.first);
3065  end_it = _displaced_mortar_residual_functors.end();
3066  }
3067 
3068  mooseAssert(
3069  it != end_it,
3070  "No ComputeMortarFunctor exists for the specified master-slave boundary pair, master "
3071  << mortar_interface.first.first << " and slave " << mortar_interface.first.second);
3072  it->second();
3073  }
3074 }
3075 
3076 void
3078 {
3079  // go over mortar constraints
3080  const auto & mortar_interfaces = _fe_problem.getMortarInterfaces(displaced);
3081 
3082  std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor<JACOBIAN>>::iterator
3083  it,
3084  end_it;
3085 
3086  for (const auto & mortar_interface : mortar_interfaces)
3087  {
3088  if (!displaced)
3089  {
3090  it = _undisplaced_mortar_jacobian_functors.find(mortar_interface.first);
3092  }
3093  else
3094  {
3095  it = _displaced_mortar_jacobian_functors.find(mortar_interface.first);
3096  end_it = _displaced_mortar_jacobian_functors.end();
3097  }
3098 
3099  mooseAssert(
3100  it != end_it,
3101  "No ComputeMortarFunctor exists for the specified master-slave boundary pair, master "
3102  << mortar_interface.first.first << " and slave " << mortar_interface.first.second);
3103  it->second();
3104  }
3105 }
virtual void setSolutionUDotDotOld(const NumericVector< Number > &u_dotdot_old)
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
Definition: SystemBase.C:805
virtual void residualSetup(THREAD_ID tid=0) const
NumericVector< Number > & getResidualTimeVector()
Return a numeric vector that is associated with the time tag.
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...
std::vector< dof_id_type > _var_all_dof_indices
TagID _Re_time_tag
Tag for time contribution residual.
virtual void addCachedResidualDirectly(NumericVector< Number > &residual, THREAD_ID tid)
Allows for all the residual contributions that are currently cached to be added directly into the vec...
ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:737
std::shared_ptr< MooseObject > create(const std::string &obj_name, const std::string &name, InputParameters parameters, THREAD_ID tid=0, bool print_deprecated=true)
Build an object (must be registered) - THIS METHOD IS DEPRECATED (Use create<T>()) ...
Definition: Factory.C:87
virtual void setSolutionUDotDot(const NumericVector< Number > &udotdot)
Set transient term used by residual and Jacobian evaluation.
virtual std::map< TagName, TagID > & getVectorTags()
Return all vector tags, where a tag is represented by a map from name to ID.
Definition: SubProblem.h:107
void reinitIncrementAtNodeForDampers(THREAD_ID tid, const std::set< MooseVariable *> &damped_vars)
Compute the incremental change in variables at nodes for dampers.
virtual void disassociateVectorFromTag(NumericVector< Number > &vec, TagID tag)
Associate a vector for a given tag.
Definition: SystemBase.C:784
Base class for deriving general dampers.
Definition: GeneralDamper.h:26
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor< ComputeStage::RESIDUAL > > _undisplaced_mortar_residual_functors
Functors for computing residuals from undisplaced mortar constraints.
MoosePreconditioner const * getPreconditioner() const
virtual const char * what() const
Get out the error message.
void setupDampers()
Setup damping stuff (called before we actually start)
virtual void addJacobianOffDiagScalar(unsigned int ivar, THREAD_ID tid=0)
void zeroVectorForResidual(const std::string &vector_name)
virtual void zeroTaggedVectors(const std::set< TagID > &tags)
Zero all vectors for given tags.
Definition: SystemBase.C:586
Base class for split-based preconditioners.
Definition: Split.h:29
bool hasActiveBlockObjects(THREAD_ID tid=0) const
bool hasVector(const std::string &tag_name) const
Check if the named vector exists in the system.
Definition: SystemBase.C:701
virtual bool uDotDotOldRequested()
Get boolean flag to check whether old solution second time derivative needs to be stored...
virtual void timestepSetup()
virtual TagID addVectorTag(TagName tag_name)
Create a Tag.
Definition: SubProblem.C:66
void setConstraintSlaveValues(NumericVector< Number > &solution, bool displaced)
Sets the value of constrained variables in the solution vector.
bool _debugging_residuals
true if debugging residuals
NumericVector< Number > * _Re_non_time
residual vector for non-time contributions
unsigned int TagID
Definition: MooseTypes.h:162
virtual NonlinearSolver< Number > * nonlinearSolver()=0
Real computeDamping(const NumericVector< Number > &solution, const NumericVector< Number > &update)
Compute damping.
virtual void setPreviousNewtonSolution(const NumericVector< Number > &soln)
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:2267
MooseObjectTagWarehouse< DGKernelBase > _dg_kernels
virtual void init() override
Initialize the system.
void reinitIncrementAtQpsForDampers(THREAD_ID tid, const std::set< MooseVariable *> &damped_vars)
Compute the incremental change in variables at QPs for dampers.
void addImplicitGeometricCouplingEntriesToJacobian(bool add=true)
If called with true this will add entries into the jacobian to link together degrees of freedom that ...
virtual NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const ParallelType type)
Adds a solution length vector to the system.
Definition: SystemBase.C:542
virtual void setSolution(const NumericVector< Number > &soln)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
Data structure used to hold penetration information.
const std::vector< std::shared_ptr< NodalConstraint > > & getActiveNodalConstraints() const
Access methods for active objects.
void addDiracKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
Adds a Dirac kernel.
bool _has_nodalbc_diag_save_in
If there is a nodal BC having diag_save_in.
virtual bool uDotRequested()
Get boolean flag to check whether solution time derivative needs to be stored.
Finds the nearest node to each node in boundary1 to each node in boundary2 and the other way around...
virtual void getDiracElements(std::set< const Elem *> &elems) override
Fills "elems" with the elements that should be looped over for Dirac Kernels.
NumericVector< Number > * _u_dotdot_old
old solution vector for u^dotdot
void checkKernelCoverage(const std::set< SubdomainID > &mesh_subdomains) const
virtual bool uDotDotRequested()
Get boolean flag to check whether solution second time derivative needs to be stored.
void addDotVectors()
Add u_dot, u_dotdot, u_dot_old and u_dotdot_old vectors if requested by the time integrator.
void getNodeDofs(dof_id_type node_id, std::vector< dof_id_type > &dofs)
Base class for all Constraint types.
Definition: Constraint.h:39
std::set< TagID > _nl_vector_tags
Vector tags to temporarily store all tags associated with the current system.
size_t _max_var_n_dofs_per_node
Maximum number of dofs for any one variable on any one node.
Definition: SystemBase.h:779
Assembly & assembly(THREAD_ID tid) override
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
virtual void associateVectorToTag(NumericVector< Number > &vec, TagID tag)
Associate a vector for a given tag.
Definition: SystemBase.C:773
MooseObjectTagWarehouse< NodalBCBase > _nodal_bcs
void addBoundaryCondition(const std::string &bc_name, const std::string &name, InputParameters parameters)
Adds a boundary condition.
virtual void setException(const std::string &message)
Set an exception.
const std::vector< std::shared_ptr< NodeFaceConstraint > > & getActiveNodeFaceConstraints(BoundaryID boundary_id, bool displaced) const
void setDecomposition(const std::vector< std::string > &decomposition)
If called with a single string, it is used as the name of a the top-level decomposition split...
bool needSubdomainMaterialOnSide(SubdomainID subdomain_id, THREAD_ID tid) const
Indicates whether this system needs material properties on internal sides.
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor< ComputeStage::JACOBIAN > > _undisplaced_mortar_jacobian_functors
Functors for computing jacobians from undisplaced mortar constraints.
Factory & _factory
Definition: SystemBase.h:738
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:472
const ElementPairInfo & getElemPairInfo(std::pair< const Elem *, const Elem *> elem_pair) const
MooseObjectWarehouseBase< Split > _splits
Decomposition splits.
NumericVector< Number > & _serialized_solution
Serialized version of the solution vector.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
bool _has_nodalbc_save_in
If there is a nodal BC having save_in.
bool hasActiveNodalConstraints() const
Deterimine if active objects exist.
MooseObjectWarehouse< NodalDamper > _nodal_dampers
Nodal Dampers for each thread.
void computeTimeDerivatives()
Computes the time derivative vector.
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 void setSolutionUDotOld(const NumericVector< Number > &u_dot_old)
bool hasDiagSaveIn() const
Weather or not the nonlinear system has diagonal Jacobian save-ins.
static PetscErrorCode Vec Mat Mat pc
void updateActive(THREAD_ID tid=0) override
Update the various active lists.
void findImplicitGeometricCouplingEntries(GeometricSearchData &geom_search_data, std::map< dof_id_type, std::vector< dof_id_type >> &graph)
Finds the implicit sparsity graph between geometrically related dofs.
void setPCSide(MooseEnum pcs)
Moose::MooseKSPNormType _ksp_norm
KSP norm type.
virtual void updateActive(THREAD_ID tid=0) override
Update the active status of Kernels.
Base class for a system (of equations)
Definition: SystemBase.h:92
void computeJacobian(SparseMatrix< Number > &jacobian, const std::set< TagID > &tags)
Associate jacobian to systemMatrixTag, and then form a matrix for all the tags.
bool computingNonlinearResid() const
Check whether residual being evaulated is non-linear.
Definition: SubProblem.h:577
std::map< dof_id_type, PenetrationInfo * > & _penetration_info
Data structure of nodes and their associated penetration information.
size_t _max_var_n_dofs_per_elem
Maximum number of dofs for any one variable on any one element.
Definition: SystemBase.h:776
bool hasDampers()
Whether or not this system has dampers.
ConstNodeRange * getLocalNodeRange()
Definition: MooseMesh.C:774
virtual void addExtraVectors() override
Method called during initialSetup to add extra system vector if they are required by the simulation...
void computeResidualTags(const std::set< TagID > &tags)
Form multiple tag-associated residual vectors for all the given tags.
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:436
TagID _Ke_system_tag
Tag for system contribution Jacobian.
void addConstraint(const std::string &c_name, const std::string &name, InputParameters parameters)
Adds a Constraint.
std::vector< std::pair< MooseVariableFEBase *, MooseVariableFEBase * > > & couplingEntries(THREAD_ID tid)
virtual TagID systemMatrixTag() override
Return the Matrix Tag ID for System.
bool _have_decomposition
Whether or not the system can be decomposed into splits.
virtual void addResidualScalar(THREAD_ID tid=0)
Scope guard for starting and stopping Floating Point Exception Trapping.
Serves as a base class for DGKernel and ADDGKernel.
Definition: DGKernelBase.h:46
void constraintResiduals(NumericVector< Number > &residual, bool displaced)
Add residual contributions from Constraints.
Base class for MOOSE preconditioners.
virtual void reinitNodesNeighbor(const std::vector< dof_id_type > &nodes, THREAD_ID tid) override
virtual GeometricSearchData & geomSearchData() override
MooseObjectTagWarehouse< KernelBase > _ad_jacobian_kernels
Base class for automatic differentiation nodal BCs that (pre)set the solution vector entries...
Specialization for filling multiple "small" preconditioning matrices simulatenously.
virtual void update()
Update the system (doing libMesh magic)
Definition: SystemBase.C:1020
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 unsigned int nVariables() const
Get the number of variables in this system.
Definition: SystemBase.C:692
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.
std::vector< std::string > split(const std::string &str, const std::string &delimiter)
Python like split function for strings.
Definition: MooseUtils.C:736
Grab all the local dof indices for the variables passed in, in the system passed in.
virtual const std::string & name() const
Definition: SystemBase.C:1088
void onTimestepBegin()
Called at the beginning of the time step.
void closeTaggedMatrices(const std::set< TagID > &tags)
Close all matrices associated the tags.
Definition: SystemBase.C:827
virtual TagID residualVectorTag() override
virtual void deactiveAllMatrixTags()
Make matrices inactive.
Definition: SystemBase.C:886
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.
NumericVector< Number > & _residual_copy
Copy of the residual vector.
void setPredictor(std::shared_ptr< Predictor > predictor)
virtual void addCachedJacobian(THREAD_ID tid) override
bool _doing_dg
true if DG is active (optimization reasons)
virtual DofMap & dofMap()
Gets writeable reference to the dof map.
Definition: SystemBase.C:932
Use whatever we have in PETSc.
Definition: MooseTypes.h:571
void addDamper(const std::string &damper_name, const std::string &name, InputParameters parameters)
Adds a damper.
bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const
Indicated whether this system needs material properties on boundaries.
FEProblemBase & _fe_problem
virtual void restoreSolutions() override
Restore current solutions (call after your solve failed)
std::vector< dof_id_type > _slave_nodes
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...
virtual std::shared_ptr< DisplacedProblem > getDisplacedProblem()
PerfID _compute_residual_tags_timer
Timers.
std::shared_ptr< T > getActiveObject(const std::string &name, THREAD_ID tid=0) const
This is the common base class for the two main kernel types implemented in MOOSE, EigenKernel and Ker...
Definition: KernelBase.h:44
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:2567
void computeDiracContributions(const std::set< TagID > &tags, bool is_jacobian)
std::vector< Point > _elem1_constraint_q_point
void updateActive(THREAD_ID tid)
Update active objects of Warehouses owned by NonlinearSystemBase.
const ElementPairList & getElemPairs() const
void setVariableGlobalDoFs(const std::string &var_name)
set all the global dof indices for a nonlinear variable
const Elem * _elem
Use whatever we have in PETSc.
Definition: MooseTypes.h:583
virtual void cacheJacobian(THREAD_ID tid) override
boundary_id_type BoundaryID
virtual void addCachedResidual(THREAD_ID tid) override
virtual void serializeSolution()
void computeJacobianTags(const std::set< TagID > &tags)
Computes multiple (tag associated) Jacobian matricese.
std::shared_ptr< MoosePreconditioner > _preconditioner
Preconditioner.
virtual void timestepSetup(THREAD_ID tid=0) const
virtual void predictorCleanup(NumericVector< Number > &ghosted_solution)
Perform cleanup tasks after application of predictor to solution vector.
DofMap & dof_map
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:112
An inteface for the _console for outputting to the Console object.
const std::vector< numeric_index_type > & n_nz
std::shared_ptr< Split > getSplit(const std::string &name)
Retrieves a split by name.
virtual void associateMatrixToTag(SparseMatrix< Number > &matrix, TagID tag)
associate a matirx to a tag
Definition: SystemBase.C:835
MooseObjectTagWarehouse< ScalarKernel > _scalar_kernels
void mortarJacobianConstraints(bool displaced)
Do mortar constraint jacobian computation.
PerfLog perf_log
Perflog to be used by applications.
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.
virtual bool uDotOldRequested()
Get boolean flag to check whether old solution time derivative needs to be stored.
void addImplicitGeometricCouplingEntries(GeometricSearchData &geom_search_data)
Adds entries to the Jacobian in the correct positions for couplings coming from dofs being coupled th...
void needsPreviousNewtonIteration(bool state)
Set a flag that indicated that user required values for the previous Newton iterate.
NumericVector< Number > & solution() override
MooseObjectTagWarehouse< KernelBase > _kernels
virtual void initialSetup()
std::map< unsigned int, std::shared_ptr< ElementPairLocator > > _element_pair_locators
MooseObjectWarehouse< T > & getMatrixTagObjectWarehouse(TagID tag_id, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has the given matrix tag...
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
Base class for deriving nodal dampers.
Definition: NodalDamper.h:32
Base class for nodal BCs that (pre)set the solution vector entries.
Definition: PresetNodalBC.h:22
unsigned int _side_num
bool areCoupled(unsigned int ivar, unsigned int jvar)
void mooseDeprecated(Args &&... args)
Emit a deprecated code/feature message with the given stringified, concatenated args.
Definition: MooseError.h:236
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.
void addSplit(const std::string &split_name, const std::string &name, InputParameters parameters)
Adds a split.
subdomain_id_type SubdomainID
virtual TagID timeVectorTag() override
Ideally, we should not need this API.
virtual void activeAllMatrixTags()
Make all exsiting matrices ative.
Definition: SystemBase.C:897
virtual unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:926
bool _need_serialized_solution
Whether or not a copy of the residual needs to be made.
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
virtual void closeTaggedVectors(const std::set< TagID > &tags)
Close all vectors for given tags.
Definition: SystemBase.C:576
virtual void reinitNodeFace(const Node *node, BoundaryID bnd_id, THREAD_ID tid) override
Base class for deriving any boundary condition that works at nodes.
Definition: NodalBCBase.h:32
NonlinearSystemBase(FEProblemBase &problem, System &sys, const std::string &name)
void computeResidualInternal(const std::set< TagID > &tags)
Compute the residual for a given tag.
Interface for objects that needs transient capabilities.
virtual std::map< TagName, TagID > & getMatrixTags()
Return all matrix tags in the sytem, where a tag is represented by a map from name to ID...
Definition: SubProblem.h:162
virtual unsigned int numVectorTags() const
The total number of tags.
Definition: SubProblem.h:122
virtual void clearDiracInfo() override
Gets called before Dirac Kernels are asked to add the points they are supposed to be evaluated in...
virtual void reinitNodes(const std::vector< dof_id_type > &nodes, THREAD_ID tid) override
virtual void jacobianSetup(THREAD_ID tid=0) const
ConstraintWarehouse _constraints
Constraints storage object.
void turnOffJacobian()
Turn off the Jacobian (must be called before equation system initialization)
MooseObjectTagWarehouse< NodalKernel > _nodal_kernels
NodalKernels for each thread.
std::map< std::pair< unsigned int, unsigned int >, PenetrationLocator * > _penetration_locators
const std::map< BoundaryID, std::vector< std::shared_ptr< T > > > & getActiveBoundaryObjects(THREAD_ID tid=0) const
virtual void addNodalKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
Adds a NodalKernel.
NumericVector< Number > & solution() override
void cacheJacobianBlock(DenseMatrix< Number > &jac_block, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, Real scaling_factor, TagID tag=0)
Definition: Assembly.C:3003
Base class for deriving element dampers.
Definition: ElementDamper.h:33
bool _add_implicit_geometric_coupling_entries_to_jacobian
Whether or not to add implicit geometric couplings to the Jacobian for FDP.
virtual void disassociateMatrixFromTag(SparseMatrix< Number > &matrix, TagID tag)
disassociate a matirx from a tag
Definition: SystemBase.C:847
virtual unsigned int numMatrixTags() const
The total number of tags.
Definition: SubProblem.h:157
Base class for creating new types of boundary conditions.
MooseApp & _app
Definition: SystemBase.h:737
void addDGKernel(std::string dg_kernel_name, const std::string &name, InputParameters parameters)
Adds a DG kernel.
std::vector< VariableWarehouse > _vars
Variable warehouses (one for each thread)
Definition: SystemBase.h:745
const std::vector< std::shared_ptr< NodeElemConstraint > > & getActiveNodeElemConstraints(SubdomainID slave_id, SubdomainID master_id, bool displaced) const
Provides a way for users to bail out of the current solve.
MatType type
unsigned int _n_residual_evaluations
Total number of residual evaluations that have been performed.
virtual void addKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
Adds a kernel.
NumericVector< Number > * _u_dot_old
old solution vector for u^dot
const SubdomainID ANY_BLOCK_ID
Definition: MooseTypes.C:15
virtual void addJacobianScalar(THREAD_ID tid=0)
virtual void prepareAssembly(THREAD_ID tid) override
virtual void reinitScalars(THREAD_ID tid) override
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.
MooseObjectWarehouse< ADPresetNodalBC< RESIDUAL > > _ad_preset_nodal_bcs
Base class for time integrators.
void addScalarKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
Adds a scalar kernel.
void constraintJacobians(bool displaced)
Add jacobian contributions from Constraints.
bool hasSaveIn() const
Weather or not the nonlinear system has save-ins.
StoredRange< std::set< const Elem * >::const_iterator, const Elem * > DistElemRange
MooseMesh & _mesh
Definition: SystemBase.h:740
bool _need_residual_copy
Whether or not a copy of the residual needs to be made.
bool hasActiveObjects(THREAD_ID tid=0) const
std::map< std::pair< unsigned int, unsigned int >, NearestNodeLocator * > _nearest_node_locators
NumericVector< Number > * _solution_previous_nl
Solution vector of the previous nonlinear iterate.
MooseObjectWarehouse< ElementDamper > _element_dampers
Element Dampers for each thread.
virtual System & system() override
Get the reference to the libMesh system.
void addObject(std::shared_ptr< Constraint > object, THREAD_ID tid=0, bool recurse=true) override
Add Constraint object to the warehouse.
virtual void setCurrentSubdomainID(const Elem *elem, THREAD_ID tid) override
PetscInt n
NumericVector< Number > * _u_dotdot
solution vector for u^dotdot
virtual void setSolutionUDot(const NumericVector< Number > &udot)
Set transient term used by residual and Jacobian evaluation.
virtual SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
Definition: SystemBase.C:811
Moose::PCSideType _pc_side
Preconditioning side.
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)
virtual void cacheJacobianNeighbor(THREAD_ID tid) override
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
Definition: SystemBase.h:782
Base class for deriving dampers.
Definition: Damper.h:29
const NumericVector< Number > * _current_solution
solution vector from nonlinear solver
virtual void reinitElemPhys(const Elem *elem, const std::vector< Point > &phys_points_in_elem, THREAD_ID tid, bool suppress_displaced_init=false) override
const std::vector< VariableName > & getVariableNames() const
Definition: SystemBase.h:710
virtual void updateActive(THREAD_ID tid=0) override
Update the active status of Kernels.
bool hasActiveElemElemConstraints(const InterfaceID interface_id, bool displaced) const
Base class for deriving any boundary condition of a integrated type.
NumericVector< Number > * _residual_ghosted
ghosted form of the residual
TagID _Re_tag
Used for the residual vector from PETSc.
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.
Moose::CouplingType coupling()
void addCachedJacobianContributions()
Adds previously-cached Jacobian values via SparseMatrix::add() calls.
Definition: Assembly.C:3488
bool hasActiveNodeElemConstraints(SubdomainID slave_id, SubdomainID master_id, bool displaced) const
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
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. ...
Definition: MooseMesh.C:2153
NumericVector< Number > * _increment_vec
increment vector
InterfaceKernelBase is the base class for all InterfaceKernel type classes.
bool doingDG() const
Getter for _doing_dg.
virtual NumericVector< Number > & serializedSolution() override
Returns a reference to a serialized version of the solution vector for this subproblem.
void computeResidualTag(NumericVector< Number > &residual, TagID tag_id)
Computes residual for a given tag.
std::vector< Point > _elem2_constraint_q_point
void addInterfaceKernel(std::string interface_kernel_name, const std::string &name, InputParameters parameters)
Adds an interface kernel.
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor< ComputeStage::RESIDUAL > > _displaced_mortar_residual_functors
Functors for computing residuals from displaced mortar constraints.
void addTimeIntegrator(const std::string &type, const std::string &name, InputParameters parameters) override
Add a time integrator.
Definition: Moose.h:112
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
virtual void cacheResidual(THREAD_ID tid) override
virtual void reinitOffDiagScalars(THREAD_ID tid) override
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.
StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode * > ConstBndNodeRange
Some useful StoredRange typedefs.
Definition: MooseMesh.h:1258
std::unordered_map< std::pair< BoundaryID, BoundaryID >, ComputeMortarFunctor< ComputeStage::JACOBIAN > > _displaced_mortar_jacobian_functors
Functors for computing jacobians from displaced mortar constraints.
std::string _decomposition_split
Name of the top-level split of the decomposition.
virtual void reinitNode(const Node *node, THREAD_ID tid) override
MooseObjectTagWarehouse< InterfaceKernelBase > _interface_kernels
NumericVector< Number > & residualVector(TagID tag)
Return a residual vector that is associated with the residual tag.
NumericVector< Number > & getResidualNonTimeVector()
Return a numeric vector that is associated with the nontime tag.
Base class for creating new types of boundary conditions.
Definition: NodalKernel.h:48
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 subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
Called from assembling when we hit a new subdomain.
virtual void restoreSolutions()
Restore current solutions (call after your solve failed)
Definition: SystemBase.C:1068
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector.
Definition: SystemBase.C:751
virtual const std::set< BoundaryID > & boundaryIDs() const
Return the boundary IDs for this object.
MooseObjectTagWarehouse< IntegratedBCBase > _integrated_bcs
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 subdomainSetup(THREAD_ID tid=0) const
StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode * > * getBoundaryNodeRange()
Definition: MooseMesh.C:788
A DiracKernel is used when you need to add contributions to the residual by means of multiplying some...
Definition: DiracKernel.h:45
virtual NumericVector< Number > & residualCopy() override
const std::vector< numeric_index_type > & n_oz
virtual NumericVector< Number > & residualGhosted() override
MooseObjectTagWarehouse< DiracKernel > _dirac_kernels
Dirac Kernel storage for each thread.
void mortarResidualConstraints(bool displaced)
Do mortar constraint residual computation.
void setMooseKSPNormType(MooseEnum kspnorm)
virtual void residualEnd(THREAD_ID tid=0) const
unsigned int THREAD_ID
Definition: MooseTypes.h:161
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:690
virtual void setNeighborSubdomainID(const Elem *elem, unsigned int side, THREAD_ID tid) override
virtual void setResidual(NumericVector< Number > &residual, THREAD_ID tid) override
const std::set< SubdomainID > & meshSubdomains() const
Returns a read-only reference to the set of subdomains currently present in the Mesh.
Definition: MooseMesh.C:2316
NumericVector< Number > * _u_dot
solution vector for u^dot
void setCachedJacobianContributions()
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
Definition: Assembly.C:3452
virtual void augmentSparsity(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.
virtual void cacheResidualNeighbor(THREAD_ID tid) override
virtual void reinitNeighborPhys(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points, THREAD_ID tid) override
MooseObjectWarehouse< PresetNodalBC > _preset_nodal_bcs
void setPreconditioner(std::shared_ptr< MoosePreconditioner > pc)
Sets a preconditioner.