www.mooseframework.org
FEProblemBase.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "FEProblemBase.h"
16 #include "AuxiliarySystem.h"
18 #include "MooseEnum.h"
19 #include "Resurrector.h"
20 #include "Factory.h"
21 #include "MooseUtils.h"
22 #include "DisplacedProblem.h"
23 #include "SystemBase.h"
24 #include "MaterialData.h"
29 #include "ComputeIndicatorThread.h"
30 #include "ComputeMarkerThread.h"
33 #include "MaxQpsThread.h"
34 #include "ActionWarehouse.h"
35 #include "Conversion.h"
36 #include "Material.h"
37 #include "ConstantIC.h"
38 #include "Parser.h"
39 #include "ElementH1Error.h"
40 #include "Function.h"
41 #include "NonlinearSystem.h"
42 #include "Distribution.h"
43 #include "Sampler.h"
44 #include "PetscSupport.h"
45 #include "RandomInterface.h"
46 #include "RandomData.h"
47 #include "MooseEigenSystem.h"
48 #include "MooseParsedFunction.h"
49 #include "MeshChangedInterface.h"
51 #include "ScalarInitialCondition.h"
52 #include "ElementPostprocessor.h"
53 #include "NodalPostprocessor.h"
54 #include "SidePostprocessor.h"
56 #include "GeneralPostprocessor.h"
62 #include "Indicator.h"
63 #include "Marker.h"
64 #include "MultiApp.h"
65 #include "MultiAppTransfer.h"
66 #include "TransientMultiApp.h"
67 #include "ElementUserObject.h"
68 #include "NodalUserObject.h"
69 #include "SideUserObject.h"
70 #include "InternalSideUserObject.h"
71 #include "GeneralUserObject.h"
72 #include "InternalSideIndicator.h"
73 #include "Transfer.h"
74 #include "MultiAppTransfer.h"
75 #include "MultiMooseEnum.h"
76 #include "Predictor.h"
77 #include "Assembly.h"
78 #include "Control.h"
79 #include "XFEMInterface.h"
80 #include "ConsoleUtils.h"
81 #include "NonlocalKernel.h"
82 #include "NonlocalIntegratedBC.h"
83 #include "ShapeElementUserObject.h"
84 #include "ShapeSideUserObject.h"
85 #include "MooseVariableScalar.h"
86 
87 #include "libmesh/exodusII_io.h"
88 #include "libmesh/quadrature.h"
89 #include "libmesh/coupling_matrix.h"
90 #include "libmesh/nonlinear_solver.h"
91 
92 // Anonymous namespace for helper function
93 namespace
94 {
98 bool
99 sortMooseVariables(MooseVariable * a, MooseVariable * b)
100 {
101  return a->number() < b->number();
102 }
103 }
104 
105 Threads::spin_mutex get_function_mutex;
106 
107 template <>
110 {
112  params.addPrivateParam<MooseMesh *>("mesh");
113  params.addParam<unsigned int>("null_space_dimension", 0, "The dimension of the nullspace");
114  params.addParam<unsigned int>(
115  "transpose_null_space_dimension", 0, "The dimension of the transpose nullspace");
116  params.addParam<unsigned int>(
117  "near_null_space_dimension", 0, "The dimension of the near nullspace");
118  params.addParam<bool>("solve",
119  true,
120  "Whether or not to actually solve the Nonlinear system. "
121  "This is handy in the case that all you want to do is "
122  "execute AuxKernels, Transfers, etc. without actually "
123  "solving anything");
124  params.addParam<bool>("use_nonlinear",
125  true,
126  "Determines whether to use a Nonlinear vs a "
127  "Eigenvalue system (Automatically determined based "
128  "on executioner)");
129  params.addParam<bool>("error_on_jacobian_nonzero_reallocation",
130  false,
131  "This causes PETSc to error if it had to reallocate memory in the Jacobian "
132  "matrix due to not having enough nonzeros");
133  params.addParam<bool>("ignore_zeros_in_jacobian",
134  false,
135  "Do not explicitly store zero values in "
136  "the Jacobian matrix if true");
137  params.addParam<bool>("force_restart",
138  false,
139  "EXPERIMENTAL: If true, a sub_app may use a "
140  "restart file instead of using of using the master "
141  "backup file");
142  params.addParam<bool>("skip_additional_restart_data",
143  false,
144  "True to skip additional data in equation system for restart. It is useful "
145  "for starting a transient calculation with a steady-state solution");
146 
147  return params;
148 }
149 
151  : SubProblem(parameters),
152  Restartable(parameters, "FEProblemBase", this),
153  _mesh(*parameters.get<MooseMesh *>("mesh")),
154  _eq(_mesh),
155  _initialized(false),
156  _kernel_type(Moose::KT_ALL),
157  _solve(getParam<bool>("solve")),
158 
159  _transient(false),
160  _time(declareRestartableData<Real>("time")),
161  _time_old(declareRestartableData<Real>("time_old")),
162  _t_step(declareRecoverableData<int>("t_step")),
163  _dt(declareRestartableData<Real>("dt")),
164  _dt_old(declareRestartableData<Real>("dt_old")),
165 
166  _nl(NULL),
167  _aux(NULL),
168  _coupling(Moose::COUPLING_DIAG),
169  _distributions(/*threaded=*/false),
170  _scalar_ics(/*threaded=*/false),
171  _material_props(
172  declareRestartableDataWithContext<MaterialPropertyStorage>("material_props", &_mesh)),
173  _bnd_material_props(
174  declareRestartableDataWithContext<MaterialPropertyStorage>("bnd_material_props", &_mesh)),
175  _pps_data(*this),
176  _vpps_data(*this),
177  _general_user_objects(/*threaded=*/false),
178  _transfers(/*threaded=*/false),
179  _to_multi_app_transfers(/*threaded=*/false),
180  _from_multi_app_transfers(/*threaded=*/false),
181 #ifdef LIBMESH_ENABLE_AMR
182  _adaptivity(*this),
183  _cycles_completed(0),
184 #endif
185  _displaced_mesh(NULL),
186  _geometric_search_data(*this, _mesh),
187  _reinit_displaced_elem(false),
188  _reinit_displaced_face(false),
189  _input_file_saved(false),
190  _has_dampers(false),
191  _has_constraints(false),
192  _has_initialized_stateful(false),
193  _const_jacobian(false),
194  _has_jacobian(false),
195  _needs_old_newton_iter(false),
196  _has_nonlocal_coupling(false),
197  _calculate_jacobian_in_uo(false),
198  _kernel_coverage_check(false),
199  _material_coverage_check(false),
200  _max_qps(std::numeric_limits<unsigned int>::max()),
201  _max_shape_funcs(std::numeric_limits<unsigned int>::max()),
202  _max_scalar_order(INVALID_ORDER),
203  _has_time_integrator(false),
204  _has_exception(false),
205  _current_execute_on_flag(EXEC_NONE),
206  _error_on_jacobian_nonzero_reallocation(
207  getParam<bool>("error_on_jacobian_nonzero_reallocation")),
208  _ignore_zeros_in_jacobian(getParam<bool>("ignore_zeros_in_jacobian")),
209  _force_restart(getParam<bool>("force_restart")),
210  _skip_additional_restart_data(getParam<bool>("skip_additional_restart_data")),
211  _fail_next_linear_convergence_check(false),
212  _currently_computing_jacobian(false),
213  _started_initial_setup(false)
214 {
215 
216  _time = 0.0;
217  _time_old = 0.0;
218  _t_step = 0;
219  _dt = 0;
220  _dt_old = _dt;
221 
222  unsigned int n_threads = libMesh::n_threads();
223 
224  _real_zero.resize(n_threads, 0.);
225  _zero.resize(n_threads);
226  _grad_zero.resize(n_threads);
227  _second_zero.resize(n_threads);
228  _second_phi_zero.resize(n_threads);
229  _uo_jacobian_moose_vars.resize(n_threads);
230 
231  _material_data.resize(n_threads);
232  _bnd_material_data.resize(n_threads);
233  _neighbor_material_data.resize(n_threads);
234  for (unsigned int i = 0; i < n_threads; i++)
235  {
236  _material_data[i] = std::make_shared<MaterialData>(_material_props);
237  _bnd_material_data[i] = std::make_shared<MaterialData>(_bnd_material_props);
238  _neighbor_material_data[i] = std::make_shared<MaterialData>(_bnd_material_props);
239  }
240 
241  _active_elemental_moose_variables.resize(n_threads);
242 
243  _block_mat_side_cache.resize(n_threads);
244  _bnd_mat_side_cache.resize(n_threads);
245 
246  _resurrector = libmesh_make_unique<Resurrector>(*this);
247 
248  _eq.parameters.set<FEProblemBase *>("_fe_problem_base") = this;
249 }
250 
251 void
253 {
254  unsigned int n_threads = libMesh::n_threads();
255 
256  _assembly.resize(n_threads);
257  for (unsigned int i = 0; i < n_threads; ++i)
258  _assembly[i] = new Assembly(nl, i);
259 }
260 
261 void
263 {
264  unsigned int n_threads = libMesh::n_threads();
265  for (unsigned int i = 0; i < n_threads; i++)
266  {
267  delete _assembly[i];
268  }
269 }
270 
271 void
273 {
274  unsigned int dimNullSpace = parameters.get<unsigned int>("null_space_dimension");
275  unsigned int dimTransposeNullSpace =
276  parameters.get<unsigned int>("transpose_null_space_dimension");
277  unsigned int dimNearNullSpace = parameters.get<unsigned int>("near_null_space_dimension");
278  for (unsigned int i = 0; i < dimNullSpace; ++i)
279  {
280  std::ostringstream oss;
281  oss << "_" << i;
282  // do not project, since this will be recomputed, but make it ghosted, since the near nullspace
283  // builder might march over all nodes
284  nl.addVector("NullSpace" + oss.str(), false, GHOSTED);
285  }
286  _subspace_dim["NullSpace"] = dimNullSpace;
287  for (unsigned int i = 0; i < dimTransposeNullSpace; ++i)
288  {
289  std::ostringstream oss;
290  oss << "_" << i;
291  // do not project, since this will be recomputed, but make it ghosted, since the near nullspace
292  // builder might march over all nodes
293  nl.addVector("TransposeNullSpace" + oss.str(), false, GHOSTED);
294  }
295  _subspace_dim["TransposeNullSpace"] = dimTransposeNullSpace;
296  for (unsigned int i = 0; i < dimNearNullSpace; ++i)
297  {
298  std::ostringstream oss;
299  oss << "_" << i;
300  // do not project, since this will be recomputed, but make it ghosted, since the near-nullspace
301  // builder might march over all semilocal nodes
302  nl.addVector("NearNullSpace" + oss.str(), false, GHOSTED);
303  }
304  _subspace_dim["NearNullSpace"] = dimNearNullSpace;
305 }
306 
308 {
309  // Flush the Console stream, the underlying call to Console::mooseConsole
310  // relies on a call to Output::checkInterval that has references to
311  // _time, etc. If it is not flushed here memory problems arise if you have
312  // an unflushed stream and start destructing things.
313  _console << std::flush;
314 
315  unsigned int n_threads = libMesh::n_threads();
316  for (unsigned int i = 0; i < n_threads; i++)
317  {
318  _zero[i].release();
319  _grad_zero[i].release();
320  _second_zero[i].release();
321  _second_phi_zero[i].release();
322  }
323 }
324 
327 {
328  std::map<SubdomainID, Moose::CoordinateSystemType>::iterator it = _coord_sys.find(sid);
329  if (it != _coord_sys.end())
330  return (*it).second;
331  else
332  mooseError("Requested subdomain ", sid, " does not exist.");
333 }
334 
335 void
336 FEProblemBase::setCoordSystem(const std::vector<SubdomainName> & blocks,
337  const MultiMooseEnum & coord_sys)
338 {
339  const std::set<SubdomainID> & subdomains = _mesh.meshSubdomains();
340  if (blocks.size() == 0)
341  {
342  // no blocks specified -> assume the whole domain
343  Moose::CoordinateSystemType coord_type = Moose::COORD_XYZ; // all is going to be XYZ by default
344  if (coord_sys.size() == 0)
345  ; // relax, do nothing
346  else if (coord_sys.size() == 1)
347  coord_type = Moose::stringToEnum<Moose::CoordinateSystemType>(
348  coord_sys[0]); // one system specified, the whole domain is going to have that system
349  else
350  mooseError("Multiple coordinate systems specified, but no blocks given.");
351 
352  for (const auto & sbd : subdomains)
353  _coord_sys[sbd] = coord_type;
354  }
355  else
356  {
357  // user specified 'blocks' but not coordinate systems
358  if (coord_sys.size() == 0)
359  {
360  // set all blocks to cartesian coordinate system
361  for (const auto & block : blocks)
362  {
363  SubdomainID sid = _mesh.getSubdomainID(block);
365  }
366  }
367  else if (coord_sys.size() == 1)
368  {
369  // set all blocks to the coordinate system specified by `coord_sys[0]`
370  Moose::CoordinateSystemType coord_type =
371  Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[0]);
372  for (const auto & block : blocks)
373  {
374  SubdomainID sid = _mesh.getSubdomainID(block);
375  _coord_sys[sid] = coord_type;
376  }
377  }
378  else
379  {
380  if (blocks.size() != coord_sys.size())
381  mooseError("Number of blocks and coordinate systems does not match.");
382 
383  for (unsigned int i = 0; i < blocks.size(); i++)
384  {
385  SubdomainID sid = _mesh.getSubdomainID(blocks[i]);
386  Moose::CoordinateSystemType coord_type =
387  Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[i]);
388  _coord_sys[sid] = coord_type;
389  }
390 
391  for (const auto & sid : subdomains)
392  if (_coord_sys.find(sid) == _coord_sys.end())
393  mooseError("Subdomain '" + Moose::stringify(sid) +
394  "' does not have a coordinate system specified.");
395  }
396  }
397 }
398 
399 void
401 {
402  _rz_coord_axis = rz_coord_axis;
403 }
404 
405 void
407 {
408  _nl->addExtraVectors();
409  _aux->addExtraVectors();
410 }
411 
412 void
414 {
415  Moose::perf_log.push("initialSetup()", "Setup");
416 
417  // set state flag indicating that we are in or beyond initialSetup.
418  // This can be used to throw errors in methods that _must_ be called at construction time.
419  _started_initial_setup = true;
420 
421  addExtraVectors();
422 
423  // Perform output related setups
425 
426  // Flush all output to _console that occur during construction and initialization of objects
428 
430  {
431  _resurrector->setRestartFile(_app.getRecoverFileBase());
432  if (_app.getRecoverFileSuffix() == "cpa")
433  _resurrector->setRestartSuffix("xda");
434  }
435 
437  _resurrector->restartFromFile();
438  else
439  {
440  ExodusII_IO * reader = _mesh.exReader();
441 
442  if (reader != NULL)
443  {
444  _nl->copyVars(*reader);
445  _aux->copyVars(*reader);
446  }
447  }
448 
449  // Build Refinement and Coarsening maps for stateful material projections if necessary
450  if (_adaptivity.isOn() &&
452  {
453  Moose::perf_log.push("mesh.buildRefinementAndCoarseningMaps()", "Setup");
455  Moose::perf_log.pop("mesh.buildRefinementAndCoarseningMaps()", "Setup");
456  }
457 
458  if (!_app.isRecovering())
459  {
466  {
467  if (!_app.isUltimateMaster())
468  mooseError(
469  "Doing extra refinements when restarting is NOT supported for sub-apps of a MultiApp");
470 
471  Moose::perf_log.push("Uniformly Refine Mesh", "Setup");
473  Moose::perf_log.pop("Uniformly Refine Mesh", "Setup");
474  }
475  }
476 
477  // Do this just in case things have been done to the mesh
479  _mesh.meshChanged();
480  if (_displaced_problem)
482 
483  unsigned int n_threads = libMesh::n_threads();
484 
485  // UserObject initialSetup
486  std::set<std::string> depend_objects_ic = _ics.getDependObjects();
487  std::set<std::string> depend_objects_aux = _aux->getDependObjects();
488 
489  _general_user_objects.updateDependObjects(depend_objects_ic, depend_objects_aux);
492 
493  for (THREAD_ID tid = 0; tid < n_threads; tid++)
494  {
495  _nodal_user_objects.updateDependObjects(depend_objects_ic, depend_objects_aux, tid);
497 
498  _elemental_user_objects.updateDependObjects(depend_objects_ic, depend_objects_aux, tid);
500 
501  _side_user_objects.updateDependObjects(depend_objects_ic, depend_objects_aux, tid);
503 
504  _internal_side_user_objects.updateDependObjects(depend_objects_ic, depend_objects_aux, tid);
506  }
507 
508  // check if jacobian calculation is done in userobject
509  for (THREAD_ID tid = 0; tid < n_threads; ++tid)
511  // Check whether nonlocal couling is required or not
515 
516  // Call the initialSetup methods for functions
517  for (THREAD_ID tid = 0; tid < n_threads; tid++)
518  {
520  tid); // initialize scalars so they are properly sized for use as input into ParsedFunctions
522  }
523 
524  if (!_app.isRecovering())
525  {
527 
528  for (THREAD_ID tid = 0; tid < n_threads; tid++)
529  _ics.initialSetup(tid);
530  _scalar_ics.sort();
531  projectSolution();
532  }
533 
534  // Materials
536  {
537  for (THREAD_ID tid = 0; tid < n_threads; tid++)
538  {
539  // Sort the Material objects, these will be actually computed by MOOSE in reinit methods.
540  _materials.sort(tid);
541 
542  // Call initialSetup on both Material and Material objects
544  }
545 
546  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
553  _assembly);
561  cmt(elem_range, true);
562 
565  }
566 
567  for (THREAD_ID tid = 0; tid < n_threads; tid++)
568  {
571  _markers.sort(tid);
572  _markers.initialSetup(tid);
573  }
574 
575 #ifdef LIBMESH_ENABLE_AMR
576 
577  if (!_app.isRecovering())
578  {
579  Moose::perf_log.push("initial adaptivity", "Setup");
580 
581  unsigned int n = adaptivity().getInitialSteps();
582  if (n && !_app.isUltimateMaster() && _app.isRestarting())
583  mooseError("Cannot perform initial adaptivity during restart on sub-apps of a MultiApp!");
584 
586  Moose::perf_log.pop("initial adaptivity", "Setup");
587  }
588 
589 #endif // LIBMESH_ENABLE_AMR
590 
591  if (!_app.isRecovering() && !_app.isRestarting())
592  {
593  // During initial setup the solution is copied to solution_old and solution_older
595  }
596 
597  if (!_app.isRecovering())
598  {
599  if (haveXFEM() && updateMeshXFEM())
600  _console << "XFEM updated mesh on initializaton" << std::endl;
601  }
602 
603  // Call initialSetup on the nonlinear system
604  _nl->initialSetup();
605 
606  // Auxilary variable initialSetup calls
607  _aux->initialSetup();
608 
609  _nl->setSolution(*(_nl->system().current_local_solution.get()));
610 
611  Moose::perf_log.push("Initial updateGeomSearch()", "Setup");
612  // Update the nearest node searches (has to be called after the problem is all set up)
613  // We do this here because this sets up the Element's DoFs to ghost
615  Moose::perf_log.pop("Initial updateGeomSearch()", "Setup");
616 
617  Moose::perf_log.push("Initial updateActiveSemiLocalNodeRange()", "Setup");
619  if (_displaced_mesh)
621  Moose::perf_log.pop("Initial updateActiveSemiLocalNodeRange()", "Setup");
622 
623  Moose::perf_log.push("reinit() after updateGeomSearch()", "Setup");
624  // Possibly reinit one more time to get ghosting correct
626  Moose::perf_log.pop("reinit() after updateGeomSearch()", "Setup");
627 
628  if (_displaced_mesh)
629  _displaced_problem->updateMesh();
630 
631  Moose::perf_log.push("Initial updateGeomSearch()", "Setup");
632  updateGeomSearch(); // Call all of the rest of the geometric searches
633  Moose::perf_log.pop("Initial updateGeomSearch()", "Setup");
634 
635  // Random interface objects
636  for (const auto & it : _random_data_objects)
637  it.second->updateSeeds(EXEC_INITIAL);
638 
639  if (_app.isRestarting() || _app.isRecovering())
640  {
641  if (_app.hasCachedBackup()) // This happens when this app is a sub-app and has been given a
642  // Backup
644  else
645  _resurrector->restartRestartableData();
646 
647  // We may have just clobbered initial conditions that were explicitly set
648  // In a _restart_ scenario it is completely valid to specify new initial conditions
649  // for some of the variables which should override what's coming from the restart file
650  if (!_app.isRecovering())
651  {
652  for (THREAD_ID tid = 0; tid < n_threads; tid++)
653  _ics.initialSetup(tid);
654  _scalar_ics.sort();
655  projectSolution();
656  }
657  }
658 
659  // HUGE NOTE: MultiApp initialSetup() MUST... I repeat MUST be _after_ restartable data has been
660  // restored
661 
662  // Call initialSetup on the MultiApps
664  {
665  _console << COLOR_CYAN << "Initializing MultiApps" << COLOR_DEFAULT << std::endl;
667  _console << COLOR_CYAN << "Finished Initializing MultiApps" << COLOR_DEFAULT << std::endl;
668  }
669 
670  // Call initialSetup on the transfers
674 
675  if (!_app.isRecovering())
676  {
678 
679  Moose::perf_log.push("execTransfers()", "Setup");
681  Moose::perf_log.pop("execTransfers()", "Setup");
682 
683  Moose::perf_log.push("execMultiApps()", "Setup");
685  if (!converged)
686  mooseError("failed to converge initial MultiApp");
687 
688  // We'll backup the Multiapp here
690  Moose::perf_log.pop("execMultiApps()", "Setup");
691 
692  // TODO: user object evaluation could fail.
694 
695  _aux->compute(EXEC_INITIAL);
696 
697  // The only user objects that should be computed here are the initial UOs
699 
701  }
702 
703  // Here we will initialize the stateful properties once more since they may have been updated
704  // during initialSetup by calls to computeProperties.
705  //
706  // It's really bad that we don't allow this during restart. It means that we can't add new
707  // stateful materials
708  // during restart. This is only happening because this _has_ to be below initial userobject
709  // execution.
710  // Otherwise this could be done up above... _before_ restoring restartable data... which would
711  // allow you to have
712  // this happen during restart. I honestly have no idea why this has to happen after initial user
713  // object computation.
714  // THAT is something we should fix... so I've opened this ticket: #5804
715  if (!_app.isRecovering() && !_app.isRestarting() &&
717  {
718  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
725  _assembly);
726  Threads::parallel_reduce(elem_range, cmt);
727  }
728 
729  // Control Logic
731 
732  // Scalar variables need to reinited for the initial conditions to be available for output
733  for (unsigned int tid = 0; tid < n_threads; tid++)
734  reinitScalars(tid);
735 
736  if (_displaced_mesh)
737  _displaced_problem->syncSolutions();
738 
739  // Writes all calls to _console from initialSetup() methods
741 
742  Moose::perf_log.pop("initialSetup()", "Setup");
743 
745  {
747  for (THREAD_ID tid = 0; tid < n_threads; ++tid)
748  _assembly[tid]->initNonlocalCoupling();
749  }
750 }
751 
752 void
754 {
755  unsigned int n_threads = libMesh::n_threads();
756 
757  for (THREAD_ID tid = 0; tid < n_threads; tid++)
758  {
761  }
762 
763  _aux->timestepSetup();
764  _nl->timestepSetup();
765 
766  // Random interface objects
767  for (const auto & it : _random_data_objects)
768  it.second->updateSeeds(EXEC_TIMESTEP_BEGIN);
769 
770  for (THREAD_ID tid = 0; tid < n_threads; tid++)
771  {
774  _markers.timestepSetup(tid);
775 
776  // Timestep setup of all UserObjects
781  }
783 
784  // Timestep setup of output objects
786 
789  _has_nonlocal_coupling = true;
790 }
791 
792 unsigned int
794 {
795  if (_max_qps == std::numeric_limits<unsigned int>::max())
796  mooseError("Max QPS uninitialized");
797  return _max_qps;
798 }
799 
800 unsigned int
802 {
803  if (_max_shape_funcs == std::numeric_limits<unsigned int>::max())
804  mooseError("Max shape functions uninitialized");
805  return _max_shape_funcs;
806 }
807 
808 Order
810 {
811  return _max_scalar_order;
812 }
813 
814 void
816 {
817  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
818  {
819  const KernelWarehouse & all_kernels = _nl->getKernelWarehouse();
820  const auto & kernels = all_kernels.getObjects(tid);
821  for (const auto & kernel : kernels)
822  {
823  std::shared_ptr<NonlocalKernel> nonlocal_kernel =
824  std::dynamic_pointer_cast<NonlocalKernel>(kernel);
825  if (nonlocal_kernel)
826  {
829  _nonlocal_kernels.addObject(kernel, tid);
830  }
831  }
832  const MooseObjectWarehouse<IntegratedBC> & all_integrated_bcs = _nl->getIntegratedBCWarehouse();
833  const auto & integrated_bcs = all_integrated_bcs.getObjects(tid);
834  for (const auto & integrated_bc : integrated_bcs)
835  {
836  std::shared_ptr<NonlocalIntegratedBC> nonlocal_integrated_bc =
837  std::dynamic_pointer_cast<NonlocalIntegratedBC>(integrated_bc);
838  if (nonlocal_integrated_bc)
839  {
842  _nonlocal_integrated_bcs.addObject(integrated_bc, tid);
843  }
844  }
845  }
846 }
847 
848 void
850 {
851  std::set<MooseVariable *> uo_jacobian_moose_vars;
852  const auto & e_objects = _elemental_user_objects.getActiveObjects(tid);
853  for (const auto & uo : e_objects)
854  {
855  std::shared_ptr<ShapeElementUserObject> shape_element_uo =
856  std::dynamic_pointer_cast<ShapeElementUserObject>(uo);
857  if (shape_element_uo)
858  {
859  _calculate_jacobian_in_uo = shape_element_uo->computeJacobianFlag();
860  const std::set<MooseVariable *> & mv_deps = shape_element_uo->jacobianMooseVariables();
861  uo_jacobian_moose_vars.insert(mv_deps.begin(), mv_deps.end());
862  }
863  }
864  const auto & s_objects = _side_user_objects.getActiveObjects(tid);
865  for (const auto & uo : s_objects)
866  {
867  std::shared_ptr<ShapeSideUserObject> shape_side_uo =
868  std::dynamic_pointer_cast<ShapeSideUserObject>(uo);
869  if (shape_side_uo)
870  {
871  _calculate_jacobian_in_uo = shape_side_uo->computeJacobianFlag();
872  const std::set<MooseVariable *> & mv_deps = shape_side_uo->jacobianMooseVariables();
873  uo_jacobian_moose_vars.insert(mv_deps.begin(), mv_deps.end());
874  }
875  }
876  _uo_jacobian_moose_vars[tid].assign(uo_jacobian_moose_vars.begin(), uo_jacobian_moose_vars.end());
877  std::sort(
878  _uo_jacobian_moose_vars[tid].begin(), _uo_jacobian_moose_vars[tid].end(), sortMooseVariables);
879 }
880 
881 void
882 FEProblemBase::setVariableAllDoFMap(const std::vector<MooseVariable *> moose_vars)
883 {
884  for (unsigned int i = 0; i < moose_vars.size(); ++i)
885  {
886  VariableName var_name = moose_vars[i]->name();
887  _nl->setVariableGlobalDoFs(var_name);
888  _var_dof_map[var_name] = _nl->getVariableGlobalDoFs();
889  }
890 }
891 
892 void
893 FEProblemBase::prepare(const Elem * elem, THREAD_ID tid)
894 {
895  _assembly[tid]->reinit(elem);
896 
897  _nl->prepare(tid);
898  _aux->prepare(tid);
899  _assembly[tid]->prepare();
901  _assembly[tid]->prepareNonlocal();
902 
904  {
905  _displaced_problem->prepare(_displaced_mesh->elemPtr(elem->id()), tid);
907  _displaced_problem->prepareNonlocal(tid);
908  }
909 }
910 
911 void
912 FEProblemBase::prepareFace(const Elem * elem, THREAD_ID tid)
913 {
914  _nl->prepareFace(tid, true);
915  _aux->prepareFace(tid, false);
916 
918  _displaced_problem->prepareFace(_displaced_mesh->elemPtr(elem->id()), tid);
919 }
920 
921 void
922 FEProblemBase::prepare(const Elem * elem,
923  unsigned int ivar,
924  unsigned int jvar,
925  const std::vector<dof_id_type> & dof_indices,
926  THREAD_ID tid)
927 {
928  _assembly[tid]->reinit(elem);
929 
930  _nl->prepare(tid);
931  _aux->prepare(tid);
932  _assembly[tid]->prepareBlock(ivar, jvar, dof_indices);
934  if (_nonlocal_cm(ivar, jvar) != 0)
935  {
936  MooseVariable & jv = _nl->getVariable(tid, jvar);
937  _assembly[tid]->prepareBlockNonlocal(ivar, jvar, dof_indices, jv.allDofIndices());
938  }
939 
941  {
942  _displaced_problem->prepare(_displaced_mesh->elemPtr(elem->id()), ivar, jvar, dof_indices, tid);
944  if (_nonlocal_cm(ivar, jvar) != 0)
945  {
946  MooseVariable & jv = _nl->getVariable(tid, jvar);
947  _displaced_problem->prepareBlockNonlocal(ivar, jvar, dof_indices, jv.allDofIndices(), tid);
948  }
949  }
950 }
951 
952 void
954 {
955  SubdomainID did = elem->subdomain_id();
956  _assembly[tid]->setCurrentSubdomainID(did);
958  _displaced_problem->assembly(tid).setCurrentSubdomainID(did);
959 }
960 
961 void
962 FEProblemBase::setNeighborSubdomainID(const Elem * elem, unsigned int side, THREAD_ID tid)
963 {
964  SubdomainID did = elem->neighbor_ptr(side)->subdomain_id();
965  _assembly[tid]->setCurrentNeighborSubdomainID(did);
967  _displaced_problem->assembly(tid).setCurrentNeighborSubdomainID(did);
968 }
969 
970 void
972 {
973  SubdomainID did = elem->subdomain_id();
974  _assembly[tid]->setCurrentNeighborSubdomainID(did);
976  _displaced_problem->assembly(tid).setCurrentNeighborSubdomainID(did);
977 }
978 
979 void
981 {
982  _assembly[tid]->prepare();
984  _assembly[tid]->prepareNonlocal();
985 
987  {
988  _displaced_problem->prepareAssembly(tid);
990  _displaced_problem->prepareNonlocal(tid);
991  }
992 }
993 
994 NumericVector<Number> &
996 {
997  return _nl->residualVector(type);
998 }
999 
1000 void
1002 {
1003  if (_nl->hasResidualVector(Moose::KT_TIME))
1005  if (_nl->hasResidualVector(Moose::KT_NONTIME))
1007 
1008  if (_displaced_problem)
1009  _displaced_problem->addResidual(tid);
1010 }
1011 
1012 void
1014 {
1015  if (_nl->hasResidualVector(Moose::KT_TIME))
1017  if (_nl->hasResidualVector(Moose::KT_NONTIME))
1019  if (_displaced_problem)
1020  _displaced_problem->addResidualNeighbor(tid);
1021 }
1022 
1023 void
1025 {
1026  if (_nl->hasResidualVector(Moose::KT_TIME))
1028  if (_nl->hasResidualVector(Moose::KT_NONTIME))
1030 }
1031 
1032 void
1034 {
1035  _assembly[tid]->cacheResidual();
1036  if (_displaced_problem)
1037  _displaced_problem->cacheResidual(tid);
1038 }
1039 
1040 void
1042 {
1043  _assembly[tid]->cacheResidualNeighbor();
1044  if (_displaced_problem)
1045  _displaced_problem->cacheResidualNeighbor(tid);
1046 }
1047 
1048 void
1050 {
1051  _assembly[tid]->addCachedResiduals();
1052 
1053  if (_displaced_problem)
1054  _displaced_problem->addCachedResidual(tid);
1055 }
1056 
1057 void
1058 FEProblemBase::addCachedResidualDirectly(NumericVector<Number> & residual, THREAD_ID tid)
1059 {
1060  if (_nl->hasResidualVector(Moose::KT_TIME))
1061  _assembly[tid]->addCachedResidual(residual, Moose::KT_TIME);
1062  if (_nl->hasResidualVector(Moose::KT_NONTIME))
1064 
1065  if (_displaced_problem)
1066  _displaced_problem->addCachedResidualDirectly(residual, tid);
1067 }
1068 
1069 void
1070 FEProblemBase::setResidual(NumericVector<Number> & residual, THREAD_ID tid)
1071 {
1072  _assembly[tid]->setResidual(residual);
1073  if (_displaced_problem)
1074  _displaced_problem->setResidual(residual, tid);
1075 }
1076 
1077 void
1078 FEProblemBase::setResidualNeighbor(NumericVector<Number> & residual, THREAD_ID tid)
1079 {
1080  _assembly[tid]->setResidualNeighbor(residual);
1081  if (_displaced_problem)
1082  _displaced_problem->setResidualNeighbor(residual, tid);
1083 }
1084 
1085 void
1086 FEProblemBase::addJacobian(SparseMatrix<Number> & jacobian, THREAD_ID tid)
1087 {
1088  _assembly[tid]->addJacobian(jacobian);
1090  _assembly[tid]->addJacobianNonlocal(jacobian);
1091  if (_displaced_problem)
1092  {
1093  _displaced_problem->addJacobian(jacobian, tid);
1095  _displaced_problem->addJacobianNonlocal(jacobian, tid);
1096  }
1097 }
1098 
1099 void
1100 FEProblemBase::addJacobianNeighbor(SparseMatrix<Number> & jacobian, THREAD_ID tid)
1101 {
1102  _assembly[tid]->addJacobianNeighbor(jacobian);
1103  if (_displaced_problem)
1104  _displaced_problem->addJacobianNeighbor(jacobian, tid);
1105 }
1106 
1107 void
1108 FEProblemBase::addJacobianScalar(SparseMatrix<Number> & jacobian, THREAD_ID tid /* = 0*/)
1109 {
1110  _assembly[tid]->addJacobianScalar(jacobian);
1111 }
1112 
1113 void
1114 FEProblemBase::addJacobianOffDiagScalar(SparseMatrix<Number> & jacobian,
1115  unsigned int ivar,
1116  THREAD_ID tid /* = 0*/)
1117 {
1118  _assembly[tid]->addJacobianOffDiagScalar(jacobian, ivar);
1119 }
1120 
1121 void
1123 {
1124  _assembly[tid]->cacheJacobian();
1126  _assembly[tid]->cacheJacobianNonlocal();
1127  if (_displaced_problem)
1128  {
1129  _displaced_problem->cacheJacobian(tid);
1131  _displaced_problem->cacheJacobianNonlocal(tid);
1132  }
1133 }
1134 
1135 void
1137 {
1138  _assembly[tid]->cacheJacobianNeighbor();
1139  if (_displaced_problem)
1140  _displaced_problem->cacheJacobianNeighbor(tid);
1141 }
1142 
1143 void
1144 FEProblemBase::addCachedJacobian(SparseMatrix<Number> & jacobian, THREAD_ID tid)
1145 {
1146  _assembly[tid]->addCachedJacobian(jacobian);
1147  if (_displaced_problem)
1148  _displaced_problem->addCachedJacobian(jacobian, tid);
1149 }
1150 
1151 void
1152 FEProblemBase::addJacobianBlock(SparseMatrix<Number> & jacobian,
1153  unsigned int ivar,
1154  unsigned int jvar,
1155  const DofMap & dof_map,
1156  std::vector<dof_id_type> & dof_indices,
1157  THREAD_ID tid)
1158 {
1159  _assembly[tid]->addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices);
1161  if (_nonlocal_cm(ivar, jvar) != 0)
1162  {
1163  MooseVariable & jv = _nl->getVariable(tid, jvar);
1164  _assembly[tid]->addJacobianBlockNonlocal(
1165  jacobian, ivar, jvar, dof_map, dof_indices, jv.allDofIndices());
1166  }
1167 
1168  if (_displaced_problem)
1169  {
1170  _displaced_problem->addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices, tid);
1172  if (_nonlocal_cm(ivar, jvar) != 0)
1173  {
1174  MooseVariable & jv = _nl->getVariable(tid, jvar);
1175  _displaced_problem->addJacobianBlockNonlocal(
1176  jacobian, ivar, jvar, dof_map, dof_indices, jv.allDofIndices(), tid);
1177  }
1178  }
1179 }
1180 
1181 void
1182 FEProblemBase::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
1183  unsigned int ivar,
1184  unsigned int jvar,
1185  const DofMap & dof_map,
1186  std::vector<dof_id_type> & dof_indices,
1187  std::vector<dof_id_type> & neighbor_dof_indices,
1188  THREAD_ID tid)
1189 {
1190  _assembly[tid]->addJacobianNeighbor(
1191  jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices);
1192  if (_displaced_problem)
1193  _displaced_problem->addJacobianNeighbor(
1194  jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices, tid);
1195 }
1196 
1197 void
1199 {
1200  _assembly[tid]->copyShapes(var);
1201 }
1202 
1203 void
1205 {
1206  _assembly[tid]->copyFaceShapes(var);
1207 }
1208 
1209 void
1211 {
1212  _assembly[tid]->copyNeighborShapes(var);
1213 }
1214 
1215 void
1216 FEProblemBase::addGhostedElem(dof_id_type elem_id)
1217 {
1218  if (_mesh.elemPtr(elem_id)->processor_id() != processor_id())
1219  _ghosted_elems.insert(elem_id);
1220 }
1221 
1222 void
1224 {
1225  _mesh.addGhostedBoundary(boundary_id);
1226  if (_displaced_problem)
1227  _displaced_mesh->addGhostedBoundary(boundary_id);
1228 }
1229 
1230 void
1232 {
1234 
1235  if (_displaced_problem)
1237 }
1238 
1239 void
1240 FEProblemBase::sizeZeroes(unsigned int /*size*/, THREAD_ID /*tid*/)
1241 {
1242  mooseDoOnce(mooseWarning(
1243  "This function is deprecated and no longer performs any function. Please do not call it."));
1244 }
1245 
1246 bool
1247 FEProblemBase::reinitDirac(const Elem * elem, THREAD_ID tid)
1248 {
1249  std::vector<Point> & points = _dirac_kernel_info.getPoints()[elem].first;
1250 
1251  unsigned int n_points = points.size();
1252 
1253  if (n_points)
1254  {
1255  if (n_points > _max_qps)
1256  {
1257  _max_qps = n_points;
1258 
1263  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
1264  {
1265  _zero[tid].resize(getMaxQps(), 0);
1266  _grad_zero[tid].resize(getMaxQps(), RealGradient(0.));
1267  _second_zero[tid].resize(getMaxQps(), RealTensor(0.));
1268  _second_phi_zero[tid].resize(
1269  getMaxQps(), std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));
1270  }
1271  }
1272 
1273  _assembly[tid]->reinitAtPhysical(elem, points);
1274 
1275  _nl->prepare(tid);
1276  _aux->prepare(tid);
1277 
1278  reinitElem(elem, tid);
1279  }
1280 
1281  _assembly[tid]->prepare();
1283  _assembly[tid]->prepareNonlocal();
1284 
1285  bool have_points = n_points > 0;
1286  if (_displaced_problem != NULL && (_reinit_displaced_elem))
1287  {
1288  have_points |= _displaced_problem->reinitDirac(_displaced_mesh->elemPtr(elem->id()), tid);
1290  _displaced_problem->prepareNonlocal(tid);
1291  }
1292 
1293  return have_points;
1294 }
1295 
1296 void
1297 FEProblemBase::reinitElem(const Elem * elem, THREAD_ID tid)
1298 {
1299  _nl->reinitElem(elem, tid);
1300  _aux->reinitElem(elem, tid);
1301 
1303  _displaced_problem->reinitElem(_displaced_mesh->elemPtr(elem->id()), tid);
1304 }
1305 
1306 void
1308  std::vector<Point> phys_points_in_elem,
1309  THREAD_ID tid)
1310 {
1311  _assembly[tid]->reinitAtPhysical(elem, phys_points_in_elem);
1312 
1313  _nl->prepare(tid);
1314  _aux->prepare(tid);
1315 
1316  reinitElem(elem, tid);
1317  _assembly[tid]->prepare();
1319  _assembly[tid]->prepareNonlocal();
1320 
1321  if (_displaced_problem != NULL && (_reinit_displaced_elem))
1322  {
1323  _displaced_problem->reinitElemPhys(
1324  _displaced_mesh->elemPtr(elem->id()), phys_points_in_elem, tid);
1326  _displaced_problem->prepareNonlocal(tid);
1327  }
1328 }
1329 
1330 void
1332  unsigned int side,
1333  BoundaryID bnd_id,
1334  THREAD_ID tid)
1335 {
1336  _assembly[tid]->reinit(elem, side);
1337 
1338  _nl->reinitElemFace(elem, side, bnd_id, tid);
1339  _aux->reinitElemFace(elem, side, bnd_id, tid);
1340 
1342  _displaced_problem->reinitElemFace(_displaced_mesh->elemPtr(elem->id()), side, bnd_id, tid);
1343 }
1344 
1345 void
1346 FEProblemBase::reinitNode(const Node * node, THREAD_ID tid)
1347 {
1348  _assembly[tid]->reinit(node);
1349 
1351  _displaced_problem->reinitNode(&_displaced_mesh->nodeRef(node->id()), tid);
1352 
1353  _nl->reinitNode(node, tid);
1354  _aux->reinitNode(node, tid);
1355 }
1356 
1357 void
1358 FEProblemBase::reinitNodeFace(const Node * node, BoundaryID bnd_id, THREAD_ID tid)
1359 {
1360  _assembly[tid]->reinit(node);
1361 
1363  _displaced_problem->reinitNodeFace(&_displaced_mesh->nodeRef(node->id()), bnd_id, tid);
1364 
1365  _nl->reinitNodeFace(node, bnd_id, tid);
1366  _aux->reinitNodeFace(node, bnd_id, tid);
1367 }
1368 
1369 void
1370 FEProblemBase::reinitNodes(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
1371 {
1373  _displaced_problem->reinitNodes(nodes, tid);
1374 
1375  _nl->reinitNodes(nodes, tid);
1376  _aux->reinitNodes(nodes, tid);
1377 }
1378 
1379 void
1380 FEProblemBase::reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
1381 {
1383  _displaced_problem->reinitNodesNeighbor(nodes, tid);
1384 
1385  _nl->reinitNodesNeighbor(nodes, tid);
1386  _aux->reinitNodesNeighbor(nodes, tid);
1387 }
1388 
1389 void
1391 {
1392  _assembly[tid]->reinitNodeNeighbor(node);
1393 
1395  _displaced_problem->reinitNodeNeighbor(&_displaced_mesh->nodeRef(node->id()), tid);
1396 
1397  _nl->reinitNodeNeighbor(node, tid);
1398  _aux->reinitNodeNeighbor(node, tid);
1399 }
1400 
1401 void
1403 {
1405  _displaced_problem->reinitScalars(tid);
1406 
1407  _nl->reinitScalars(tid);
1408  _aux->reinitScalars(tid);
1409 
1410  _assembly[tid]->prepareScalar();
1411 }
1412 
1413 void
1415 {
1416  _assembly[tid]->prepareOffDiagScalar();
1417  if (_displaced_problem != NULL)
1418  _displaced_problem->reinitOffDiagScalars(tid);
1419 }
1420 
1421 void
1422 FEProblemBase::reinitNeighbor(const Elem * elem, unsigned int side, THREAD_ID tid)
1423 {
1424  const Elem * neighbor = elem->neighbor_ptr(side);
1425  unsigned int neighbor_side = neighbor->which_neighbor_am_i(elem);
1426 
1427  _assembly[tid]->reinitElemAndNeighbor(elem, side, neighbor, neighbor_side);
1428 
1429  _nl->prepareNeighbor(tid);
1430  _aux->prepareNeighbor(tid);
1431 
1432  _assembly[tid]->prepareNeighbor();
1433 
1434  BoundaryID bnd_id = 0; // some dummy number (it is not really used for anything, right now)
1435  _nl->reinitElemFace(elem, side, bnd_id, tid);
1436  _aux->reinitElemFace(elem, side, bnd_id, tid);
1437 
1438  _nl->reinitNeighborFace(neighbor, neighbor_side, bnd_id, tid);
1439  _aux->reinitNeighborFace(neighbor, neighbor_side, bnd_id, tid);
1440 
1442  _displaced_problem->reinitNeighbor(elem, side, tid);
1443 }
1444 
1445 void
1447  unsigned int neighbor_side,
1448  const std::vector<Point> & physical_points,
1449  THREAD_ID tid)
1450 {
1451  // Reinits shape the functions at the physical points
1452  _assembly[tid]->reinitNeighborAtPhysical(neighbor, neighbor_side, physical_points);
1453 
1454  // Sets the neighbor dof indices
1455  _nl->prepareNeighbor(tid);
1456  _aux->prepareNeighbor(tid);
1457 
1458  // Resizes Re and Ke
1459  _assembly[tid]->prepareNeighbor();
1460 
1461  // Compute the values of each variable at the points
1462  _nl->reinitNeighborFace(neighbor, neighbor_side, 0, tid);
1463  _aux->reinitNeighborFace(neighbor, neighbor_side, 0, tid);
1464 
1465  // Do the same for the displaced problem
1466  if (_displaced_problem != NULL)
1467  _displaced_problem->reinitNeighborPhys(
1468  _displaced_mesh->elemPtr(neighbor->id()), neighbor_side, physical_points, tid);
1469 }
1470 
1471 void
1473  const std::vector<Point> & physical_points,
1474  THREAD_ID tid)
1475 {
1476  // Reinits shape the functions at the physical points
1477  _assembly[tid]->reinitNeighborAtPhysical(neighbor, physical_points);
1478 
1479  // Sets the neighbor dof indices
1480  _nl->prepareNeighbor(tid);
1481  _aux->prepareNeighbor(tid);
1482 
1483  // Resizes Re and Ke
1484  _assembly[tid]->prepareNeighbor();
1485 
1486  // Compute the values of each variable at the points
1487  _nl->reinitNeighbor(neighbor, tid);
1488  _aux->reinitNeighbor(neighbor, tid);
1489 
1490  // Do the same for the displaced problem
1491  if (_displaced_problem != NULL)
1492  _displaced_problem->reinitNeighborPhys(
1493  _displaced_mesh->elemPtr(neighbor->id()), physical_points, tid);
1494 }
1495 
1496 void
1497 FEProblemBase::getDiracElements(std::set<const Elem *> & elems)
1498 {
1499  // First add in the undisplaced elements
1500  elems = _dirac_kernel_info.getElements();
1501 
1502  if (_displaced_problem)
1503  {
1504  std::set<const Elem *> displaced_elements;
1505  _displaced_problem->getDiracElements(displaced_elements);
1506 
1507  { // Use the ids from the displaced elements to get the undisplaced elements
1508  // and add them to the list
1509  for (const auto & elem : displaced_elements)
1510  elems.insert(_mesh.elemPtr(elem->id()));
1511  }
1512  }
1513 }
1514 
1515 void
1517 {
1519 
1520  if (_displaced_problem)
1521  _displaced_problem->clearDiracInfo();
1522 }
1523 
1524 void
1526 {
1527  _all_materials.subdomainSetup(subdomain, tid);
1528 
1529  // Call the subdomain methods of the output system, these are not threaded so only call it once
1530  if (tid == 0)
1532 
1533  _nl->subdomainSetup(subdomain, tid);
1534 
1535  // FIXME: call displaced_problem->subdomainSetup() ?
1536  // When adding possibility with materials being evaluated on displaced mesh
1537 }
1538 
1539 void
1541 {
1542  _all_materials.neighborSubdomainSetup(subdomain, tid);
1543 }
1544 
1545 void
1547 {
1548  setInputParametersFEProblem(parameters);
1549  parameters.set<SubProblem *>("_subproblem") = this;
1550 
1551  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1552  {
1553  std::shared_ptr<Function> func = _factory.create<Function>(type, name, parameters, tid);
1554  _functions.addObject(func, tid);
1555  }
1556 }
1557 
1558 bool
1559 FEProblemBase::hasFunction(const std::string & name, THREAD_ID tid)
1560 {
1561  return _functions.hasActiveObject(name, tid);
1562 }
1563 
1564 Function &
1565 FEProblemBase::getFunction(const std::string & name, THREAD_ID tid)
1566 {
1567  // This thread lock is necessary since this method will create functions
1568  // for all threads if one is missing.
1569  Threads::spin_mutex::scoped_lock lock(get_function_mutex);
1570 
1571  if (!hasFunction(name, tid))
1572  {
1573  // If we didn't find a function, it might be a default function, attempt to construct one now
1574  std::istringstream ss(name);
1575  Real real_value;
1576 
1577  // First see if it's just a constant. If it is, build a ConstantFunction
1578  if (ss >> real_value && ss.eof())
1579  {
1580  InputParameters params = _factory.getValidParams("ConstantFunction");
1581  params.set<Real>("value") = real_value;
1582  addFunction("ConstantFunction", ss.str(), params);
1583  }
1584  else
1585  {
1587  std::string vars = "x,y,z,t,NaN,pi,e";
1588  if (fp.Parse(name, vars) == -1) // -1 for success
1589  {
1590  // It parsed ok, so build a MooseParsedFunction
1591  InputParameters params = _factory.getValidParams("ParsedFunction");
1592  params.set<std::string>("value") = name;
1593  addFunction("ParsedFunction", name, params);
1594  }
1595  }
1596 
1597  // Try once more
1598  if (!hasFunction(name, tid))
1599  mooseError("Unable to find function " + name);
1600  }
1601 
1602  return *(_functions.getActiveObject(name, tid));
1603 }
1604 
1607 {
1608  mooseDeprecated("FEProblemBase::getNonlinearSystem() is deprecated, please use "
1609  "FEProblemBase::getNonlinearSystemBase() \n");
1610 
1611  auto nl_sys = std::dynamic_pointer_cast<NonlinearSystem>(_nl);
1612 
1613  if (!nl_sys)
1614  mooseError("This is not a NonlinearSystem");
1615 
1616  return *nl_sys;
1617 }
1618 
1619 void
1621  const std::string & name,
1623 {
1624  setInputParametersFEProblem(parameters);
1625  parameters.set<std::string>("type") = type;
1626  std::shared_ptr<Distribution> dist = _factory.create<Distribution>(type, name, parameters);
1627  _distributions.addObject(dist);
1628 }
1629 
1630 Distribution &
1632 {
1633  if (!_distributions.hasActiveObject(name))
1634  mooseError("Unable to find distribution " + name);
1635 
1636  return *(_distributions.getActiveObject(name));
1637 }
1638 
1639 void
1641 {
1642  setInputParametersFEProblem(parameters);
1643  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1644  {
1645  std::shared_ptr<Sampler> dist = _factory.create<Sampler>(type, name, parameters, tid);
1646  _samplers.addObject(dist, tid);
1647  }
1648 }
1649 
1650 Sampler &
1651 FEProblemBase::getSampler(const std::string & name, THREAD_ID tid)
1652 {
1653  if (!_samplers.hasActiveObject(name, tid))
1654  mooseError("Unable to find Sampler " + name);
1655 
1656  return *(_samplers.getActiveObject(name, tid));
1657 }
1658 
1659 bool
1660 FEProblemBase::duplicateVariableCheck(const std::string & var_name,
1661  const FEType & type,
1662  bool is_aux)
1663 {
1664  SystemBase * curr_sys_ptr = _nl.get();
1665  SystemBase * other_sys_ptr = _aux.get();
1666  std::string error_prefix = "";
1667  if (is_aux)
1668  {
1669  curr_sys_ptr = _aux.get();
1670  other_sys_ptr = _nl.get();
1671  error_prefix = "Aux";
1672  }
1673 
1674  if (other_sys_ptr->hasVariable(var_name))
1675  mooseError("Cannot have an auxiliary variable and a nonlinear variable with the same name: ",
1676  var_name);
1677 
1678  if (curr_sys_ptr->hasVariable(var_name))
1679  {
1680  const Variable & var =
1681  curr_sys_ptr->system().variable(curr_sys_ptr->system().variable_number(var_name));
1682  if (var.type() != type)
1683  mooseError(error_prefix,
1684  "Variable with name '",
1685  var_name,
1686  "' already exists but is of a differing type!");
1687 
1688  return true;
1689  }
1690 
1691  return false;
1692 }
1693 
1694 void
1695 FEProblemBase::addVariable(const std::string & var_name,
1696  const FEType & type,
1697  Real scale_factor,
1698  const std::set<SubdomainID> * const active_subdomains)
1699 {
1700  if (duplicateVariableCheck(var_name, type, /* is_aux = */ false))
1701  return;
1702 
1703  _nl->addVariable(var_name, type, scale_factor, active_subdomains);
1704  if (_displaced_problem)
1705  _displaced_problem->addVariable(var_name, type, scale_factor, active_subdomains);
1706 }
1707 
1708 void
1709 FEProblemBase::addScalarVariable(const std::string & var_name,
1710  Order order,
1711  Real scale_factor,
1712  const std::set<SubdomainID> * const active_subdomains)
1713 {
1714  if (order > _max_scalar_order)
1715  _max_scalar_order = order;
1716 
1717  FEType type(order, SCALAR);
1718  if (duplicateVariableCheck(var_name, type, /* is_aux = */ false))
1719  return;
1720 
1721  _nl->addScalarVariable(var_name, order, scale_factor, active_subdomains);
1722  if (_displaced_problem)
1723  _displaced_problem->addScalarVariable(var_name, order, scale_factor, active_subdomains);
1724 }
1725 
1726 void
1727 FEProblemBase::addKernel(const std::string & kernel_name,
1728  const std::string & name,
1730 {
1731  setInputParametersFEProblem(parameters);
1732  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1733  {
1734  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1735  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1736  _reinit_displaced_elem = true;
1737  }
1738  else
1739  {
1740  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1741  {
1742  // We allow Kernels to request that they use_displaced_mesh,
1743  // but then be overridden when no displacements variables are
1744  // provided in the Mesh block. If that happened, update the value
1745  // of use_displaced_mesh appropriately for this Kernel.
1746  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1747  parameters.set<bool>("use_displaced_mesh") = false;
1748  }
1749 
1750  parameters.set<SubProblem *>("_subproblem") = this;
1751  parameters.set<SystemBase *>("_sys") = _nl.get();
1752  }
1753 
1754  _nl->addKernel(kernel_name, name, parameters);
1755 }
1756 
1757 void
1758 FEProblemBase::addNodalKernel(const std::string & kernel_name,
1759  const std::string & name,
1761 {
1762  setInputParametersFEProblem(parameters);
1763  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1764  {
1765  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1766  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1767  _reinit_displaced_elem = true;
1768  }
1769  else
1770  {
1771  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1772  {
1773  // We allow NodalKernels to request that they use_displaced_mesh,
1774  // but then be overridden when no displacements variables are
1775  // provided in the Mesh block. If that happened, update the value
1776  // of use_displaced_mesh appropriately for this NodalKernel.
1777  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1778  parameters.set<bool>("use_displaced_mesh") = false;
1779  }
1780 
1781  parameters.set<SubProblem *>("_subproblem") = this;
1782  parameters.set<SystemBase *>("_sys") = _nl.get();
1783  }
1784  _nl->addNodalKernel(kernel_name, name, parameters);
1785 }
1786 
1787 void
1788 FEProblemBase::addScalarKernel(const std::string & kernel_name,
1789  const std::string & name,
1791 {
1792  setInputParametersFEProblem(parameters);
1793  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1794  {
1795  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1796  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1797  }
1798  else
1799  {
1800  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1801  {
1802  // We allow ScalarKernels to request that they use_displaced_mesh,
1803  // but then be overridden when no displacements variables are
1804  // provided in the Mesh block. If that happened, update the value
1805  // of use_displaced_mesh appropriately for this ScalarKernel.
1806  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1807  parameters.set<bool>("use_displaced_mesh") = false;
1808  }
1809 
1810  parameters.set<SubProblem *>("_subproblem") = this;
1811  parameters.set<SystemBase *>("_sys") = _nl.get();
1812  }
1813  _nl->addScalarKernel(kernel_name, name, parameters);
1814 }
1815 
1816 void
1817 FEProblemBase::addBoundaryCondition(const std::string & bc_name,
1818  const std::string & name,
1820 {
1821  setInputParametersFEProblem(parameters);
1822  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1823  {
1824  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1825  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1826  _reinit_displaced_face = true;
1827  }
1828  else
1829  {
1830  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1831  {
1832  // We allow Materials to request that they use_displaced_mesh,
1833  // but then be overridden when no displacements variables are
1834  // provided in the Mesh block. If that happened, update the value
1835  // of use_displaced_mesh appropriately for this Material.
1836  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1837  parameters.set<bool>("use_displaced_mesh") = false;
1838  }
1839 
1840  parameters.set<SubProblem *>("_subproblem") = this;
1841  parameters.set<SystemBase *>("_sys") = _nl.get();
1842  }
1843  _nl->addBoundaryCondition(bc_name, name, parameters);
1844 }
1845 
1846 void
1847 FEProblemBase::addConstraint(const std::string & c_name,
1848  const std::string & name,
1850 {
1851  _has_constraints = true;
1852 
1853  setInputParametersFEProblem(parameters);
1854  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1855  {
1856  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1857  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1858  _reinit_displaced_face = true;
1859  }
1860  else
1861  {
1862  // It might _want_ to use a displaced mesh... but we're not so set it to false
1863  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1864  parameters.set<bool>("use_displaced_mesh") = false;
1865 
1866  parameters.set<SubProblem *>("_subproblem") = this;
1867  parameters.set<SystemBase *>("_sys") = _nl.get();
1868  }
1869  _nl->addConstraint(c_name, name, parameters);
1870 }
1871 
1872 void
1873 FEProblemBase::addAuxVariable(const std::string & var_name,
1874  const FEType & type,
1875  const std::set<SubdomainID> * const active_subdomains)
1876 {
1877  if (duplicateVariableCheck(var_name, type, /* is_aux = */ true))
1878  return;
1879 
1880  _aux->addVariable(var_name, type, 1.0, active_subdomains);
1881  if (_displaced_problem)
1882  _displaced_problem->addAuxVariable(var_name, type, active_subdomains);
1883 }
1884 
1885 void
1886 FEProblemBase::addAuxScalarVariable(const std::string & var_name,
1887  Order order,
1888  Real scale_factor,
1889  const std::set<SubdomainID> * const active_subdomains)
1890 {
1891  if (order > _max_scalar_order)
1892  _max_scalar_order = order;
1893 
1894  FEType type(order, SCALAR);
1895  if (duplicateVariableCheck(var_name, type, /* is_aux = */ true))
1896  return;
1897 
1898  _aux->addScalarVariable(var_name, order, scale_factor, active_subdomains);
1899  if (_displaced_problem)
1900  _displaced_problem->addAuxScalarVariable(var_name, order, scale_factor, active_subdomains);
1901 }
1902 
1903 void
1904 FEProblemBase::addAuxKernel(const std::string & kernel_name,
1905  const std::string & name,
1907 {
1908  setInputParametersFEProblem(parameters);
1909  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1910  {
1911  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1912  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
1913  parameters.set<SystemBase *>("_nl_sys") = &_displaced_problem->nlSys();
1914  if (!parameters.get<std::vector<BoundaryName>>("boundary").empty())
1915  _reinit_displaced_face = true;
1916  else
1917  _reinit_displaced_elem = true;
1918  }
1919  else
1920  {
1921  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1922  {
1923  // We allow AuxKernels to request that they use_displaced_mesh,
1924  // but then be overridden when no displacements variables are
1925  // provided in the Mesh block. If that happened, update the value
1926  // of use_displaced_mesh appropriately for this AuxKernel.
1927  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1928  parameters.set<bool>("use_displaced_mesh") = false;
1929  }
1930 
1931  parameters.set<SubProblem *>("_subproblem") = this;
1932  parameters.set<SystemBase *>("_sys") = _aux.get();
1933  parameters.set<SystemBase *>("_nl_sys") = _nl.get();
1934  }
1935  _aux->addKernel(kernel_name, name, parameters);
1936 }
1937 
1938 void
1939 FEProblemBase::addAuxScalarKernel(const std::string & kernel_name,
1940  const std::string & name,
1942 {
1943  setInputParametersFEProblem(parameters);
1944  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1945  {
1946  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1947  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
1948  }
1949  else
1950  {
1951  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1952  {
1953  // We allow AuxScalarKernels to request that they use_displaced_mesh,
1954  // but then be overridden when no displacements variables are
1955  // provided in the Mesh block. If that happened, update the value
1956  // of use_displaced_mesh appropriately for this AuxScalarKernel.
1957  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1958  parameters.set<bool>("use_displaced_mesh") = false;
1959  }
1960 
1961  parameters.set<SubProblem *>("_subproblem") = this;
1962  parameters.set<SystemBase *>("_sys") = _aux.get();
1963  }
1964  _aux->addScalarKernel(kernel_name, name, parameters);
1965 }
1966 
1967 void
1968 FEProblemBase::addDiracKernel(const std::string & kernel_name,
1969  const std::string & name,
1971 {
1972  setInputParametersFEProblem(parameters);
1973  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
1974  {
1975  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1976  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1977  _reinit_displaced_elem = true;
1978  }
1979  else
1980  {
1981  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
1982  {
1983  // We allow DiracKernels to request that they use_displaced_mesh,
1984  // but then be overridden when no displacements variables are
1985  // provided in the Mesh block. If that happened, update the value
1986  // of use_displaced_mesh appropriately for this DiracKernel.
1987  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1988  parameters.set<bool>("use_displaced_mesh") = false;
1989  }
1990 
1991  parameters.set<SubProblem *>("_subproblem") = this;
1992  parameters.set<SystemBase *>("_sys") = _nl.get();
1993  }
1994  _nl->addDiracKernel(kernel_name, name, parameters);
1995 }
1996 
1997 // DGKernels ////
1998 
1999 void
2000 FEProblemBase::addDGKernel(const std::string & dg_kernel_name,
2001  const std::string & name,
2003 {
2004  setInputParametersFEProblem(parameters);
2005  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
2006  {
2007  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2008  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2009  _reinit_displaced_face = true;
2010  }
2011  else
2012  {
2013  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
2014  {
2015  // We allow DGKernels to request that they use_displaced_mesh,
2016  // but then be overridden when no displacements variables are
2017  // provided in the Mesh block. If that happened, update the value
2018  // of use_displaced_mesh appropriately for this DGKernel.
2019  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2020  parameters.set<bool>("use_displaced_mesh") = false;
2021  }
2022 
2023  parameters.set<SubProblem *>("_subproblem") = this;
2024  parameters.set<SystemBase *>("_sys") = _nl.get();
2025  }
2026  _nl->addDGKernel(dg_kernel_name, name, parameters);
2027 }
2028 
2029 // InterfaceKernels ////
2030 
2031 void
2032 FEProblemBase::addInterfaceKernel(const std::string & interface_kernel_name,
2033  const std::string & name,
2035 {
2036  setInputParametersFEProblem(parameters);
2037  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
2038  {
2039  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2040  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2041  _reinit_displaced_face = true;
2042  }
2043  else
2044  {
2045  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
2046  {
2047  // We allow InterfaceKernels to request that they use_displaced_mesh,
2048  // but then be overridden when no displacements variables are
2049  // provided in the Mesh block. If that happened, update the value
2050  // of use_displaced_mesh appropriately for this InterfaceKernel.
2051  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2052  parameters.set<bool>("use_displaced_mesh") = false;
2053  }
2054 
2055  parameters.set<SubProblem *>("_subproblem") = this;
2056  parameters.set<SystemBase *>("_sys") = _nl.get();
2057  }
2058  _nl->addInterfaceKernel(interface_kernel_name, name, parameters);
2059 }
2060 
2061 void
2062 FEProblemBase::addInitialCondition(const std::string & ic_name,
2063  const std::string & name,
2065 {
2066 
2067  // before we start to mess with the initial condition, we need to check parameters for errors.
2068  parameters.checkParams(name);
2069 
2070  setInputParametersFEProblem(parameters);
2071  parameters.set<SubProblem *>("_subproblem") = this;
2072 
2073  const std::string & var_name = parameters.get<VariableName>("variable");
2074  // field IC
2075  if (hasVariable(var_name))
2076  {
2077  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2078  {
2079  MooseVariable & var = getVariable(tid, var_name);
2080  parameters.set<SystemBase *>("_sys") = &var.sys();
2081  std::shared_ptr<InitialCondition> ic =
2082  _factory.create<InitialCondition>(ic_name, name, parameters, tid);
2083  _ics.addObject(ic, tid);
2084  }
2085  }
2086 
2087  // scalar IC
2088  else if (hasScalarVariable(var_name))
2089  {
2090  MooseVariableScalar & var = getScalarVariable(0, var_name);
2091  parameters.set<SystemBase *>("_sys") = &var.sys();
2092  std::shared_ptr<ScalarInitialCondition> ic =
2094  _scalar_ics.addObject(ic);
2095  }
2096 
2097  else
2098  mooseError(
2099  "Variable '", var_name, "' requested in initial condition '", name, "' does not exist.");
2100 }
2101 
2102 void
2104 {
2105  Moose::perf_log.push("projectSolution()", "Utility");
2106 
2107  Moose::enableFPE();
2108 
2109  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
2110  ComputeInitialConditionThread cic(*this);
2111  Threads::parallel_reduce(elem_range, cic);
2112 
2113  // Need to close the solution vector here so that boundary ICs take precendence
2114  _nl->solution().close();
2115  _aux->solution().close();
2116 
2117  // now run boundary-restricted initial conditions
2118  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
2120  Threads::parallel_reduce(bnd_nodes, cbic);
2121 
2122  _nl->solution().close();
2123  _aux->solution().close();
2124 
2125  // Also, load values into the SCALAR dofs
2126  // Note: We assume that all SCALAR dofs are on the
2127  // processor with highest ID
2128  if (processor_id() == (n_processors() - 1) && _scalar_ics.hasActiveObjects())
2129  {
2130  const auto & ics = _scalar_ics.getActiveObjects();
2131  for (const auto & ic : ics)
2132  {
2133  MooseVariableScalar & var = ic->variable();
2134  var.reinit();
2135 
2136  DenseVector<Number> vals(var.order());
2137  ic->compute(vals);
2138 
2139  const unsigned int n_SCALAR_dofs = var.dofIndices().size();
2140  for (unsigned int i = 0; i < n_SCALAR_dofs; i++)
2141  {
2142  const dof_id_type global_index = var.dofIndices()[i];
2143  var.sys().solution().set(global_index, vals(i));
2144  var.setValue(i, vals(i));
2145  }
2146  }
2147  }
2148 
2149  Moose::enableFPE(false);
2150 
2151  _nl->solution().close();
2152  _nl->solution().localize(*_nl->system().current_local_solution, _nl->dofMap().get_send_list());
2153 
2154  _aux->solution().close();
2155  _aux->solution().localize(*_aux->sys().current_local_solution, _aux->dofMap().get_send_list());
2156 
2157  Moose::perf_log.pop("projectSolution()", "Utility");
2158 }
2159 
2160 std::shared_ptr<Material>
2163  THREAD_ID tid,
2164  bool no_warn)
2165 {
2166  switch (type)
2167  {
2169  name += "_neighbor";
2170  break;
2172  name += "_face";
2173  break;
2174  default:
2175  break;
2176  }
2177 
2178  std::shared_ptr<Material> material = _all_materials[type].getActiveObject(name, tid);
2179  if (!no_warn && material->getParam<bool>("compute") && type == Moose::BLOCK_MATERIAL_DATA)
2180  mooseWarning("You are retrieving a Material object (",
2181  material->name(),
2182  "), but its compute flag is set to true. This indicates that MOOSE is "
2183  "computing this property which may not be desired and produce un-expected "
2184  "results.");
2185 
2186  return material;
2187 }
2188 
2189 std::shared_ptr<MaterialData>
2191 {
2192  std::shared_ptr<MaterialData> output;
2193  switch (type)
2194  {
2196  output = _material_data[tid];
2197  break;
2199  output = _neighbor_material_data[tid];
2200  break;
2203  output = _bnd_material_data[tid];
2204  break;
2205  }
2206  return output;
2207 }
2208 
2209 void
2210 FEProblemBase::addMaterial(const std::string & mat_name,
2211  const std::string & name,
2213 {
2214  setInputParametersFEProblem(parameters);
2215  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
2216  {
2217  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2218  _reinit_displaced_elem = true;
2219  }
2220  else
2221  {
2222  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
2223  {
2224  // We allow Materials to request that they use_displaced_mesh,
2225  // but then be overridden when no displacements variables are
2226  // provided in the Mesh block. If that happened, update the value
2227  // of use_displaced_mesh appropriately for this Material.
2228  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2229  parameters.set<bool>("use_displaced_mesh") = false;
2230  }
2231 
2232  parameters.set<SubProblem *>("_subproblem") = this;
2233  }
2234 
2235  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2236  {
2237  // Create the general Block/Boundary Material object
2238  std::shared_ptr<Material> material = _factory.create<Material>(mat_name, name, parameters, tid);
2239  bool discrete = !material->getParam<bool>("compute");
2240 
2241  // If the object is boundary restricted do not create the neighbor and face objects
2242  if (material->boundaryRestricted())
2243  {
2244  _all_materials.addObject(material, tid);
2245  if (discrete)
2246  _discrete_materials.addObject(material, tid);
2247  else
2248  _materials.addObject(material, tid);
2249  }
2250 
2251  // Non-boundary restricted require face and neighbor objects
2252  else
2253  {
2254  // The name of the object being created, this is changed multiple times as objects are created
2255  // below
2256  std::string object_name;
2257 
2258  // Create a copy of the supplied parameters to the setting for "_material_data_type" isn't
2259  // used from a previous tid loop
2260  InputParameters current_parameters = parameters;
2261 
2262  // face material
2263  current_parameters.set<Moose::MaterialDataType>("_material_data_type") =
2265  object_name = name + "_face";
2266  std::shared_ptr<Material> face_material =
2267  _factory.create<Material>(mat_name, object_name, current_parameters, tid);
2268 
2269  // neighbor material
2270  current_parameters.set<Moose::MaterialDataType>("_material_data_type") =
2272  current_parameters.set<bool>("_neighbor") = true;
2273  object_name = name + "_neighbor";
2274  std::shared_ptr<Material> neighbor_material =
2275  _factory.create<Material>(mat_name, object_name, current_parameters, tid);
2276 
2277  // Store the material objects
2278  _all_materials.addObjects(material, neighbor_material, face_material, tid);
2279 
2280  if (discrete)
2281  _discrete_materials.addObjects(material, neighbor_material, face_material, tid);
2282  else
2283  _materials.addObjects(material, neighbor_material, face_material, tid);
2284 
2285  // link enabled parameter of face and neighbor materials
2286  MooseObjectParameterName name(MooseObjectName("Material", material->name()), "enabled");
2287  MooseObjectParameterName face_name(MooseObjectName("Material", face_material->name()),
2288  "enabled");
2289  MooseObjectParameterName neighbor_name(MooseObjectName("Material", neighbor_material->name()),
2290  "enabled");
2293  }
2294  }
2295 }
2296 
2297 void
2299 {
2300  std::set<MooseVariable *> needed_moose_vars;
2301  std::set<unsigned int> needed_mat_props;
2302 
2303  if (_all_materials.hasActiveBlockObjects(blk_id, tid))
2304  {
2305  _all_materials.updateVariableDependency(needed_moose_vars, tid);
2306  _all_materials.updateBlockMatPropDependency(blk_id, needed_mat_props, tid);
2307  }
2308 
2309  const std::set<BoundaryID> & ids = _mesh.getSubdomainBoundaryIds(blk_id);
2310  for (const auto & id : ids)
2311  {
2312  _materials.updateBoundaryVariableDependency(id, needed_moose_vars, tid);
2313  _materials.updateBoundaryMatPropDependency(id, needed_mat_props, tid);
2314  }
2315 
2316  const std::set<MooseVariable *> & current_active_elemental_moose_variables =
2318  needed_moose_vars.insert(current_active_elemental_moose_variables.begin(),
2319  current_active_elemental_moose_variables.end());
2320 
2321  const std::set<unsigned int> & current_active_material_properties =
2323  needed_mat_props.insert(current_active_material_properties.begin(),
2324  current_active_material_properties.end());
2325 
2326  setActiveElementalMooseVariables(needed_moose_vars, tid);
2327  setActiveMaterialProperties(needed_mat_props, tid);
2328 }
2329 
2330 void
2331 FEProblemBase::reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful)
2332 {
2333  if (hasActiveMaterialProperties(tid))
2334  {
2335  const Elem *& elem = _assembly[tid]->elem();
2336  unsigned int n_points = _assembly[tid]->qRule()->n_points();
2337  _material_data[tid]->resize(n_points);
2338 
2339  // Only swap if requested
2340  if (swap_stateful)
2341  _material_data[tid]->swap(*elem);
2342 
2343  if (_discrete_materials.hasActiveBlockObjects(blk_id, tid))
2344  _material_data[tid]->reset(_discrete_materials.getActiveBlockObjects(blk_id, tid));
2345 
2346  if (_materials.hasActiveBlockObjects(blk_id, tid))
2347  _material_data[tid]->reinit(_materials.getActiveBlockObjects(blk_id, tid));
2348  }
2349 }
2350 
2351 void
2353 {
2354  if (hasActiveMaterialProperties(tid))
2355  {
2356  const Elem *& elem = _assembly[tid]->elem();
2357  unsigned int side = _assembly[tid]->side();
2358  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2359 
2360  _bnd_material_data[tid]->resize(n_points);
2361 
2362  if (swap_stateful && !_bnd_material_data[tid]->isSwapped())
2363  _bnd_material_data[tid]->swap(*elem, side);
2364 
2365  if (_discrete_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2366  _bnd_material_data[tid]->reset(
2367  _discrete_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2368 
2369  if (_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2370  _bnd_material_data[tid]->reinit(
2371  _materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2372  }
2373 }
2374 
2375 void
2377 {
2378  if (hasActiveMaterialProperties(tid))
2379  {
2380  // NOTE: this will not work with h-adaptivity
2381  const Elem *& neighbor = _assembly[tid]->neighbor();
2382  unsigned int neighbor_side = neighbor->which_neighbor_am_i(_assembly[tid]->elem());
2383  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2384  _neighbor_material_data[tid]->resize(n_points);
2385 
2386  // Only swap if requested
2387  if (swap_stateful)
2388  _neighbor_material_data[tid]->swap(*neighbor, neighbor_side);
2389 
2390  if (_discrete_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2391  _neighbor_material_data[tid]->reset(
2392  _discrete_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2393 
2394  if (_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2395  _neighbor_material_data[tid]->reinit(
2396  _materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2397  }
2398 }
2399 
2400 void
2401 FEProblemBase::reinitMaterialsBoundary(BoundaryID boundary_id, THREAD_ID tid, bool swap_stateful)
2402 {
2403  if (hasActiveMaterialProperties(tid))
2404  {
2405  const Elem *& elem = _assembly[tid]->elem();
2406  unsigned int side = _assembly[tid]->side();
2407  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2408  _bnd_material_data[tid]->resize(n_points);
2409 
2410  if (swap_stateful && !_bnd_material_data[tid]->isSwapped())
2411  _bnd_material_data[tid]->swap(*elem, side);
2412 
2413  if (_discrete_materials.hasActiveBoundaryObjects(boundary_id, tid))
2414  _bnd_material_data[tid]->reset(
2415  _discrete_materials.getActiveBoundaryObjects(boundary_id, tid));
2416 
2417  if (_materials.hasActiveBoundaryObjects(boundary_id, tid))
2418  _bnd_material_data[tid]->reinit(_materials.getActiveBoundaryObjects(boundary_id, tid));
2419  }
2420 }
2421 
2422 void
2424 {
2425  const Elem *& elem = _assembly[tid]->elem();
2426  _material_data[tid]->swapBack(*elem);
2427 }
2428 
2429 void
2431 {
2432  const Elem *& elem = _assembly[tid]->elem();
2433  unsigned int side = _assembly[tid]->side();
2434  _bnd_material_data[tid]->swapBack(*elem, side);
2435 }
2436 
2437 void
2439 {
2440  // NOTE: this will not work with h-adaptivity
2441  const Elem *& neighbor = _assembly[tid]->neighbor();
2442  unsigned int neighbor_side = neighbor->which_neighbor_am_i(_assembly[tid]->elem());
2443  _neighbor_material_data[tid]->swapBack(*neighbor, neighbor_side);
2444 }
2445 
2450 std::shared_ptr<Postprocessor>
2451 getPostprocessorPointer(std::shared_ptr<MooseObject> mo)
2452 {
2453  {
2454  std::shared_ptr<ElementPostprocessor> intermediate =
2455  std::dynamic_pointer_cast<ElementPostprocessor>(mo);
2456  if (intermediate.get())
2457  return std::static_pointer_cast<Postprocessor>(intermediate);
2458  }
2459 
2460  {
2461  std::shared_ptr<NodalPostprocessor> intermediate =
2462  std::dynamic_pointer_cast<NodalPostprocessor>(mo);
2463  if (intermediate.get())
2464  return std::static_pointer_cast<Postprocessor>(intermediate);
2465  }
2466 
2467  {
2468  std::shared_ptr<InternalSidePostprocessor> intermediate =
2469  std::dynamic_pointer_cast<InternalSidePostprocessor>(mo);
2470  if (intermediate.get())
2471  return std::static_pointer_cast<Postprocessor>(intermediate);
2472  }
2473 
2474  {
2475  std::shared_ptr<SidePostprocessor> intermediate =
2476  std::dynamic_pointer_cast<SidePostprocessor>(mo);
2477  if (intermediate.get())
2478  return std::static_pointer_cast<Postprocessor>(intermediate);
2479  }
2480 
2481  {
2482  std::shared_ptr<GeneralPostprocessor> intermediate =
2483  std::dynamic_pointer_cast<GeneralPostprocessor>(mo);
2484  if (intermediate.get())
2485  return std::static_pointer_cast<Postprocessor>(intermediate);
2486  }
2487 
2488  return std::shared_ptr<Postprocessor>();
2489 }
2490 
2491 template <typename UO_TYPE, typename PP_TYPE>
2492 Postprocessor *
2494 {
2495  PP_TYPE * intermediate = dynamic_cast<PP_TYPE *>(uo);
2496  if (intermediate)
2497  return static_cast<Postprocessor *>(intermediate);
2498 
2499  return NULL;
2500 }
2501 
2502 void
2504 {
2505  _pps_data.init(name);
2506 }
2507 
2508 void
2510  const std::string & name,
2512 {
2513  // Check for name collision
2515  mooseError(std::string("A UserObject with the name \"") + name +
2516  "\" already exists. You may not add a Postprocessor by the same name.");
2517 
2518  addUserObject(pp_name, name, parameters);
2519  initPostprocessorData(name);
2520 }
2521 
2522 void
2524  const std::string & name,
2526 {
2527  // Check for name collision
2529  mooseError(std::string("A UserObject with the name \"") + name +
2530  "\" already exists. You may not add a VectorPostprocessor by the same name.");
2531 
2532  addUserObject(pp_name, name, parameters);
2533 }
2534 
2535 void
2536 FEProblemBase::addUserObject(std::string user_object_name,
2537  const std::string & name,
2539 {
2540  setInputParametersFEProblem(parameters);
2541  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
2542  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2543  else
2544  {
2545  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
2546  {
2547  // We allow UserObjects to request that they use_displaced_mesh,
2548  // but then be overridden when no displacements variables are
2549  // provided in the Mesh block. If that happened, update the value
2550  // of use_displaced_mesh appropriately for this UserObject.
2551  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2552  parameters.set<bool>("use_displaced_mesh") = false;
2553  }
2554 
2555  parameters.set<SubProblem *>("_subproblem") = this;
2556  }
2557 
2558  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2559  {
2560  // Create the UserObject
2561  std::shared_ptr<UserObject> user_object =
2562  _factory.create<UserObject>(user_object_name, name, parameters, tid);
2563  _all_user_objects.addObject(user_object, tid);
2564 
2565  // Attempt to create all the possible UserObject types
2566  std::shared_ptr<ElementUserObject> euo =
2567  std::dynamic_pointer_cast<ElementUserObject>(user_object);
2568  std::shared_ptr<SideUserObject> suo = std::dynamic_pointer_cast<SideUserObject>(user_object);
2569  std::shared_ptr<InternalSideUserObject> isuo =
2570  std::dynamic_pointer_cast<InternalSideUserObject>(user_object);
2571  std::shared_ptr<NodalUserObject> nuo = std::dynamic_pointer_cast<NodalUserObject>(user_object);
2572  std::shared_ptr<GeneralUserObject> guo =
2573  std::dynamic_pointer_cast<GeneralUserObject>(user_object);
2574 
2575  // Account for displaced mesh use
2576  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
2577  {
2578  if (euo || nuo)
2579  _reinit_displaced_elem = true;
2580  else if (suo)
2581  _reinit_displaced_face = true;
2582  }
2583 
2584  // Add the object to the correct warehouse
2585  if (guo)
2586  {
2588  break; // not threaded
2589  }
2590  else if (nuo)
2591  _nodal_user_objects.addObject(nuo, tid);
2592  else if (suo)
2593  _side_user_objects.addObject(suo, tid);
2594  else if (isuo)
2596  else if (euo)
2598  }
2599 }
2600 
2601 const UserObject &
2603 {
2605  return *(_all_user_objects.getActiveObject(name).get());
2606 
2607  mooseError("Unable to find user object with name '" + name + "'");
2608 }
2609 
2610 bool
2612 {
2613  return _all_user_objects.hasActiveObject(name);
2614 }
2615 
2616 bool
2618 {
2619  return _pps_data.hasPostprocessor(name);
2620 }
2621 
2623 FEProblemBase::getPostprocessorValue(const PostprocessorName & name)
2624 {
2625  return _pps_data.getPostprocessorValue(name);
2626 }
2627 
2630 {
2631  return _pps_data.getPostprocessorValueOld(name);
2632 }
2633 
2636 {
2638 }
2639 
2642 {
2643  return _vpps_data;
2644 }
2645 
2646 bool
2648 {
2649  return _vpps_data.hasVectorPostprocessor(name);
2650 }
2651 
2653 FEProblemBase::getVectorPostprocessorValue(const VectorPostprocessorName & name,
2654  const std::string & vector_name)
2655 {
2656  return _vpps_data.getVectorPostprocessorValue(name, vector_name);
2657 }
2658 
2661  const std::string & vector_name)
2662 {
2663  return _vpps_data.getVectorPostprocessorValueOld(name, vector_name);
2664 }
2665 
2668  const std::string & vector_name)
2669 {
2670  return _vpps_data.declareVector(name, vector_name);
2671 }
2672 
2673 const std::vector<std::pair<std::string, VectorPostprocessorData::VectorPostprocessorState>> &
2675 {
2676  return _vpps_data.vectors(vpp_name);
2677 }
2678 
2679 void
2681 {
2682  for (const auto & it : _multi_apps)
2683  {
2684  const auto & objects = it.second.getActiveObjects();
2685  for (const auto & obj : objects)
2686  obj->parentOutputPositionChanged();
2687  }
2688 }
2689 
2690 void
2692 {
2694  computeMarkers();
2695 }
2696 
2697 void
2699 {
2700  // Initialize indicator aux variable fields
2702  {
2703  Moose::perf_log.push("Adaptivity: computeIndicators()", "Execution");
2704 
2705  std::vector<std::string> fields;
2706 
2707  // Indicator Fields
2708  const auto & indicators = _indicators.getActiveObjects();
2709  for (const auto & indicator : indicators)
2710  fields.push_back(indicator->name());
2711 
2712  // InternalSideIndicator Fields
2713  const auto & internal_indicators = _internal_side_indicators.getActiveObjects();
2714  for (const auto & internal_indicator : internal_indicators)
2715  fields.push_back(internal_indicator->name());
2716 
2717  _aux->zeroVariables(fields);
2718 
2719  // compute Indicators
2720  ComputeIndicatorThread cit(*this);
2721  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cit);
2722  _aux->solution().close();
2723  _aux->update();
2724 
2725  ComputeIndicatorThread finalize_cit(*this, true);
2726  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), finalize_cit);
2727  _aux->solution().close();
2728  _aux->update();
2729 
2730  Moose::perf_log.pop("Adaptivity: computeIndicators()", "Execution");
2731  }
2732 }
2733 
2734 void
2736 {
2737  if (_markers.hasActiveObjects())
2738  {
2739  Moose::perf_log.push("Adaptivity: computeMarkers()", "Execution");
2740 
2741  std::vector<std::string> fields;
2742 
2743  // Marker Fields
2744  const auto & markers = _markers.getActiveObjects();
2745  for (const auto & marker : markers)
2746  fields.push_back(marker->name());
2747 
2748  _aux->zeroVariables(fields);
2749 
2751 
2752  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2753  {
2754  const auto & markers = _markers.getActiveObjects(tid);
2755  for (const auto & marker : markers)
2756  marker->markerSetup();
2757  }
2758 
2759  ComputeMarkerThread cmt(*this);
2760  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cmt);
2761 
2762  _aux->solution().close();
2763  _aux->update();
2764 
2765  Moose::perf_log.pop("Adaptivity: computeMarkers()", "Execution");
2766  }
2767 }
2768 
2769 const ExecFlagType &
2771 {
2772  return _current_execute_on_flag;
2773 }
2774 
2775 void
2777 {
2778  // Set the current flag
2779  _current_execute_on_flag = exec_type;
2780  if (exec_type == EXEC_NONLINEAR)
2782 
2783  // Samplers
2784  executeSamplers(exec_type);
2785 
2786  // Pre-aux UserObjects
2787  computeUserObjects(exec_type, Moose::PRE_AUX);
2788 
2789  // AuxKernels
2790  computeAuxiliaryKernels(exec_type);
2791 
2792  // Post-aux UserObjects
2793  computeUserObjects(exec_type, Moose::POST_AUX);
2794 
2795  // Controls
2796  executeControls(exec_type);
2797 
2798  // Return the current flag to None
2801 }
2802 
2803 void
2805 {
2806  _aux->compute(type);
2807 }
2808 
2809 void
2811 {
2812  // Get convenience reference to active warehouse
2815  const MooseObjectWarehouse<InternalSideUserObject> & internal_side =
2819 
2820  if (!elemental.hasActiveObjects() && !side.hasActiveObjects() &&
2821  !internal_side.hasActiveObjects() && !nodal.hasActiveObjects() && !general.hasActiveObjects())
2822  // Nothing to do, return early
2823  return;
2824 
2825  // Start the timer here since we have at least one active user object
2826  std::string compute_uo_tag = "computeUserObjects(" + Moose::stringify(type) + ")";
2827  Moose::perf_log.push(compute_uo_tag, "Execution");
2828 
2829  // Perform Residual/Jacobian setups
2830  switch (type)
2831  {
2832  case EXEC_LINEAR:
2833  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2834  {
2835  elemental.residualSetup(tid);
2836  side.residualSetup(tid);
2837  internal_side.residualSetup(tid);
2838  nodal.residualSetup(tid);
2839  }
2840  general.residualSetup();
2841  break;
2842 
2843  case EXEC_NONLINEAR:
2844  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2845  {
2846  elemental.jacobianSetup(tid);
2847  side.jacobianSetup(tid);
2848  internal_side.jacobianSetup(tid);
2849  nodal.jacobianSetup(tid);
2850  }
2851  general.jacobianSetup();
2852  break;
2853 
2854  default:
2855  break;
2856  }
2857 
2858  // Initialize Elemental/Side/InternalSideUserObjects
2859  initializeUserObjects<ElementUserObject>(elemental);
2860  initializeUserObjects<SideUserObject>(side);
2861  initializeUserObjects<InternalSideUserObject>(internal_side);
2862 
2863  // Execute Elemental/Side/InternalSideUserObjects
2864  if (elemental.hasActiveObjects() || side.hasActiveObjects() || internal_side.hasActiveObjects())
2865  {
2866  ComputeUserObjectsThread cppt(*this, getNonlinearSystemBase(), elemental, side, internal_side);
2867  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cppt);
2868  }
2869 
2870  // Finalize, threadJoin, and update PP values of Elemental/Side/InternalSideUserObjects
2871  finalizeUserObjects<SideUserObject>(side);
2872  finalizeUserObjects<InternalSideUserObject>(internal_side);
2873  finalizeUserObjects<ElementUserObject>(elemental);
2874 
2875  // Initialize Nodal
2876  initializeUserObjects<NodalUserObject>(nodal);
2877 
2878  // Execute NodalUserObjects
2879  if (nodal.hasActiveObjects())
2880  {
2881  ComputeNodalUserObjectsThread cnppt(*this, nodal);
2882  Threads::parallel_reduce(*_mesh.getLocalNodeRange(), cnppt);
2883  }
2884 
2885  // Finalize, threadJoin, and update PP values of Nodal
2886  finalizeUserObjects<NodalUserObject>(nodal);
2887 
2888  // Execute GeneralUserObjects
2889  if (general.hasActiveObjects())
2890  {
2891  const auto & objects = general.getActiveObjects();
2892  for (const auto & obj : objects)
2893  {
2894  obj->initialize();
2895  obj->execute();
2896  obj->finalize();
2897 
2898  std::shared_ptr<Postprocessor> pp = std::dynamic_pointer_cast<Postprocessor>(obj);
2899  if (pp)
2900  _pps_data.storeValue(obj->name(), pp->getValue());
2901  }
2902  }
2903 
2904  Moose::perf_log.pop(compute_uo_tag, "Execution");
2905 }
2906 
2907 void
2909 {
2911 
2912  auto controls_wh = _control_warehouse[exec_type];
2913  // Add all of the dependencies into the resolver and sort them
2914  for (const auto & it : controls_wh.getActiveObjects())
2915  {
2916  // Make sure an item with no dependencies comes out too!
2917  resolver.addItem(it);
2918 
2919  std::vector<std::string> & dependent_controls = it->getDependencies();
2920  for (const auto & depend_name : dependent_controls)
2921  {
2922  if (controls_wh.hasActiveObject(depend_name))
2923  {
2924  auto dep_control = controls_wh.getActiveObject(depend_name);
2925  resolver.insertDependency(it, dep_control);
2926  }
2927  else
2928  mooseError("The Control \"",
2929  depend_name,
2930  "\" was not created, did you make a "
2931  "spelling mistake or forget to include it "
2932  "in your input file?");
2933  }
2934  }
2935 
2936  const auto & ordered_controls = resolver.getSortedValues();
2937 
2938  if (!ordered_controls.empty())
2939  {
2940  Moose::perf_log.push("computeControls()", "Execution");
2941 
2942  _control_warehouse.setup(exec_type);
2943  // Run the controls in the proper order
2944  for (const auto & control : ordered_controls)
2945  control->execute();
2946 
2947  Moose::perf_log.pop("computeControls()", "Execution");
2948  }
2949 }
2950 
2951 void
2953 {
2954  // TODO: This should be done in a threaded loop, but this should be super quick so for now
2955  // do a serial loop.
2956  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2957  {
2958  const auto & objects = _samplers[exec_type].getActiveObjects(tid);
2959  if (!objects.empty())
2960  {
2961  Moose::perf_log.push("executeSamplers()", "Execution");
2962  _samplers.setup(exec_type);
2963  for (auto & sampler : objects)
2964  sampler->execute();
2965  Moose::perf_log.pop("executeSamplers()", "Execution");
2966  }
2967  }
2968 }
2969 
2970 void
2972 {
2973  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2974  {
2975  _nl->updateActive(tid);
2976  _aux->updateActive(tid);
2979  _markers.updateActive(tid);
2981  _materials.updateActive(tid);
2987  _samplers.updateActive(tid);
2988  }
2989 
2997 }
2998 
2999 void
3001 {
3002  //<< "Object " << a->name() << " -> " << b->name() << std::endl;
3003 }
3004 
3005 void
3007 {
3008  // Need to see if _any_ processor has ghosted elems or geometry objects.
3009  bool needs_reinit = !_ghosted_elems.empty();
3010  needs_reinit = needs_reinit || !_geometric_search_data._nearest_node_locators.empty();
3011  needs_reinit =
3012  needs_reinit ||
3013  (_displaced_problem && !_displaced_problem->geomSearchData()._nearest_node_locators.empty());
3014  _communicator.max(needs_reinit);
3015 
3016  if (needs_reinit)
3017  {
3018  // Call reinit to get the ghosted vectors correct now that some geometric search has been done
3019  _eq.reinit();
3020 
3021  if (_displaced_mesh)
3022  _displaced_problem->es().reinit();
3023  }
3024 }
3025 
3026 void
3027 FEProblemBase::addDamper(std::string damper_name,
3028  const std::string & name,
3030 {
3031  setInputParametersFEProblem(parameters);
3032  parameters.set<SubProblem *>("_subproblem") = this;
3033  parameters.set<SystemBase *>("_sys") = _nl.get();
3034 
3035  _has_dampers = true;
3036  _nl->addDamper(damper_name, name, parameters);
3037 }
3038 
3039 void
3041 {
3042  _nl->setupDampers();
3043 }
3044 
3045 void
3046 FEProblemBase::addIndicator(std::string indicator_name,
3047  const std::string & name,
3049 {
3050  setInputParametersFEProblem(parameters);
3051  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
3052  {
3053  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3054  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3055  _reinit_displaced_elem = true;
3056  }
3057  else
3058  {
3059  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
3060  {
3061  // We allow Indicators to request that they use_displaced_mesh,
3062  // but then be overridden when no displacements variables are
3063  // provided in the Mesh block. If that happened, update the value
3064  // of use_displaced_mesh appropriately for this Indicator.
3065  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3066  parameters.set<bool>("use_displaced_mesh") = false;
3067  }
3068 
3069  parameters.set<SubProblem *>("_subproblem") = this;
3070  parameters.set<SystemBase *>("_sys") = _aux.get();
3071  }
3072 
3073  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
3074  {
3075  std::shared_ptr<Indicator> indicator =
3076  _factory.create<Indicator>(indicator_name, name, parameters, tid);
3077 
3078  std::shared_ptr<InternalSideIndicator> isi =
3079  std::dynamic_pointer_cast<InternalSideIndicator>(indicator);
3080  if (isi)
3082  else
3083  _indicators.addObject(indicator, tid);
3084  }
3085 }
3086 
3087 void
3088 FEProblemBase::addMarker(std::string marker_name,
3089  const std::string & name,
3091 {
3092  setInputParametersFEProblem(parameters);
3093  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
3094  {
3095  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3096  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3097  _reinit_displaced_elem = true;
3098  }
3099  else
3100  {
3101  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
3102  {
3103  // We allow Markers to request that they use_displaced_mesh,
3104  // but then be overridden when no displacements variables are
3105  // provided in the Mesh block. If that happened, update the value
3106  // of use_displaced_mesh appropriately for this Marker.
3107  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3108  parameters.set<bool>("use_displaced_mesh") = false;
3109  }
3110 
3111  parameters.set<SubProblem *>("_subproblem") = this;
3112  parameters.set<SystemBase *>("_sys") = _aux.get();
3113  }
3114 
3115  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
3116  {
3117  std::shared_ptr<Marker> marker = _factory.create<Marker>(marker_name, name, parameters, tid);
3118  _markers.addObject(marker, tid);
3119  }
3120 }
3121 
3122 void
3123 FEProblemBase::addMultiApp(const std::string & multi_app_name,
3124  const std::string & name,
3126 {
3127  setInputParametersFEProblem(parameters);
3128  parameters.set<MPI_Comm>("_mpi_comm") = _communicator.get();
3129  parameters.set<std::shared_ptr<CommandLine>>("_command_line") = _app.commandLine();
3130 
3131  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
3132  {
3133  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3134  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3135  _reinit_displaced_elem = true;
3136  }
3137  else
3138  {
3139  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
3140  {
3141  // We allow MultiApps to request that they use_displaced_mesh,
3142  // but then be overridden when no displacements variables are
3143  // provided in the Mesh block. If that happened, update the value
3144  // of use_displaced_mesh appropriately for this MultiApp.
3145  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3146  parameters.set<bool>("use_displaced_mesh") = false;
3147  }
3148 
3149  parameters.set<SubProblem *>("_subproblem") = this;
3150  parameters.set<SystemBase *>("_sys") = _aux.get();
3151  }
3152 
3153  std::shared_ptr<MultiApp> multi_app = _factory.create<MultiApp>(multi_app_name, name, parameters);
3154 
3155  _multi_apps.addObject(multi_app);
3156 
3157  // Store TranseintMultiApp objects in another container, this is needed for calling computeDT
3158  std::shared_ptr<TransientMultiApp> trans_multi_app =
3159  std::dynamic_pointer_cast<TransientMultiApp>(multi_app);
3160  if (trans_multi_app)
3161  _transient_multi_apps.addObject(trans_multi_app);
3162 }
3163 
3164 bool
3165 FEProblemBase::hasMultiApp(const std::string & multi_app_name)
3166 {
3167  return _multi_apps.hasActiveObject(multi_app_name);
3168 }
3169 
3170 std::shared_ptr<MultiApp>
3171 FEProblemBase::getMultiApp(const std::string & multi_app_name)
3172 {
3173  return _multi_apps.getActiveObject(multi_app_name);
3174 }
3175 
3176 void
3178 {
3179  bool to_multiapp = direction == MultiAppTransfer::TO_MULTIAPP;
3180  std::string string_direction = to_multiapp ? " To " : " From ";
3181  const MooseObjectWarehouse<Transfer> & wh =
3183 
3184  if (wh.hasActiveObjects())
3185  {
3186  const auto & transfers = wh.getActiveObjects();
3187 
3188  _console << COLOR_CYAN << "\nStarting Transfers on " << Moose::stringify(type)
3189  << string_direction << "MultiApps" << COLOR_DEFAULT << std::endl;
3190  for (const auto & transfer : transfers)
3191  {
3192  Moose::perf_log.push(transfer->name(), "Transfers");
3193  transfer->execute();
3194  Moose::perf_log.pop(transfer->name(), "Transfers");
3195  }
3196 
3197  _console << "Waiting For Transfers To Finish" << '\n';
3198  MooseUtils::parallelBarrierNotify(_communicator);
3199 
3200  _console << COLOR_CYAN << "Transfers on " << Moose::stringify(type) << " Are Finished\n"
3201  << COLOR_DEFAULT << std::endl;
3202  }
3203  else if (_multi_apps[type].getActiveObjects().size())
3204  _console << COLOR_CYAN << "\nNo Transfers on " << Moose::stringify(type) << " To MultiApps\n"
3205  << COLOR_DEFAULT << std::endl;
3206 }
3207 
3208 std::vector<std::shared_ptr<Transfer>>
3210 {
3214  return wh.getActiveObjects();
3215 }
3216 
3217 bool
3219 {
3220  // Active MultiApps
3221  const std::vector<MooseSharedPointer<MultiApp>> & multi_apps =
3223 
3224  // Do anything that needs to be done to Apps before transfers
3225  for (const auto & multi_app : multi_apps)
3226  multi_app->preTransfer(_dt, _time);
3227 
3228  // Execute Transfers _to_ MultiApps
3230 
3231  // Execute MultiApps
3232  if (multi_apps.size())
3233  {
3234  _console << COLOR_CYAN << "\nExecuting MultiApps on " << Moose::stringify(type) << COLOR_DEFAULT
3235  << std::endl;
3236 
3237  bool success = true;
3238 
3239  for (const auto & multi_app : multi_apps)
3240  success = multi_app->solveStep(_dt, _time, auto_advance);
3241 
3242  _console << "Waiting For Other Processors To Finish" << '\n';
3243  MooseUtils::parallelBarrierNotify(_communicator);
3244 
3245  _communicator.max(success);
3246 
3247  if (!success)
3248  return false;
3249 
3250  _console << COLOR_CYAN << "Finished Executing MultiApps on " << Moose::stringify(type) << "\n"
3251  << COLOR_DEFAULT << std::endl;
3252  }
3253 
3254  // Execute Transfers _from_ MultiApps
3256 
3257  // If we made it here then everything passed
3258  return true;
3259 }
3260 
3261 void
3263 {
3264  const auto & multi_apps = _multi_apps.getActiveObjects();
3265 
3266  for (const auto & multi_app : multi_apps)
3267  // If the app has been solved, then postExecute() will have been called already too
3268  if (!multi_app->isSolved())
3269  multi_app->postExecute();
3270 }
3271 
3272 void
3274 {
3275  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3276 
3277  if (multi_apps.size())
3278  {
3279  _console << COLOR_CYAN << "\nAdvancing MultiApps" << COLOR_DEFAULT << std::endl;
3280 
3281  for (const auto & multi_app : multi_apps)
3282  multi_app->advanceStep();
3283 
3284  _console << "Waiting For Other Processors To Finish" << std::endl;
3285  MooseUtils::parallelBarrierNotify(_communicator);
3286 
3287  _console << COLOR_CYAN << "Finished Advancing MultiApps\n" << COLOR_DEFAULT << std::endl;
3288  }
3289 }
3290 
3291 void
3293 {
3294  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3295 
3296  if (multi_apps.size())
3297  {
3298  _console << COLOR_CYAN << "\nBacking Up MultiApps" << COLOR_DEFAULT << std::endl;
3299 
3300  for (const auto & multi_app : multi_apps)
3301  multi_app->backup();
3302 
3303  _console << "Waiting For Other Processors To Finish" << std::endl;
3304  MooseUtils::parallelBarrierNotify(_communicator);
3305 
3306  _console << COLOR_CYAN << "Finished Backing Up MultiApps\n" << COLOR_DEFAULT << std::endl;
3307  }
3308 }
3309 
3310 void
3312 {
3313  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3314 
3315  if (multi_apps.size())
3316  {
3317  if (force)
3318  _console << COLOR_CYAN << "\nRestoring Multiapps because of solve failure!" << COLOR_DEFAULT
3319  << std::endl;
3320  else
3321  _console << COLOR_CYAN << "\nRestoring MultiApps" << COLOR_DEFAULT << std::endl;
3322 
3323  for (const auto & multi_app : multi_apps)
3324  if (force || multi_app->needsRestoration())
3325  multi_app->restore();
3326 
3327  _console << "Waiting For Other Processors To Finish" << std::endl;
3328  MooseUtils::parallelBarrierNotify(_communicator);
3329 
3330  _console << COLOR_CYAN << "Finished Restoring MultiApps\n" << COLOR_DEFAULT << std::endl;
3331  }
3332 }
3333 
3334 Real
3336 {
3337  const auto & multi_apps = _transient_multi_apps[type].getActiveObjects();
3338 
3339  Real smallest_dt = std::numeric_limits<Real>::max();
3340 
3341  for (const auto & multi_app : multi_apps)
3342  smallest_dt = std::min(smallest_dt, multi_app->computeDT());
3343 
3344  return smallest_dt;
3345 }
3346 
3347 void
3349 {
3350  if (_transfers[type].hasActiveObjects())
3351  {
3352  const auto & transfers = _transfers[type].getActiveObjects();
3353  for (const auto & transfer : transfers)
3354  transfer->execute();
3355  }
3356 }
3357 
3358 void
3359 FEProblemBase::addTransfer(const std::string & transfer_name,
3360  const std::string & name,
3362 {
3363  setInputParametersFEProblem(parameters);
3364  if (_displaced_problem != NULL && parameters.get<bool>("use_displaced_mesh"))
3365  {
3366  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3367  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3368  _reinit_displaced_elem = true;
3369  }
3370  else
3371  {
3372  if (_displaced_problem == NULL && parameters.get<bool>("use_displaced_mesh"))
3373  {
3374  // We allow Transfers to request that they use_displaced_mesh,
3375  // but then be overridden when no displacements variables are
3376  // provided in the Mesh block. If that happened, update the value
3377  // of use_displaced_mesh appropriately for this Transfer.
3378  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3379  parameters.set<bool>("use_displaced_mesh") = false;
3380  }
3381 
3382  parameters.set<SubProblem *>("_subproblem") = this;
3383  parameters.set<SystemBase *>("_sys") = _aux.get();
3384  }
3385 
3386  // Create the Transfer objects
3387  std::shared_ptr<Transfer> transfer = _factory.create<Transfer>(transfer_name, name, parameters);
3388 
3389  // Add MultiAppTransfer object
3390  std::shared_ptr<MultiAppTransfer> multi_app_transfer =
3391  std::dynamic_pointer_cast<MultiAppTransfer>(transfer);
3392  if (multi_app_transfer)
3393  {
3394  if (multi_app_transfer->direction() == MultiAppTransfer::TO_MULTIAPP)
3395  _to_multi_app_transfers.addObject(multi_app_transfer);
3396  else
3397  _from_multi_app_transfers.addObject(multi_app_transfer);
3398  }
3399  else
3400  _transfers.addObject(transfer);
3401 }
3402 
3403 bool
3404 FEProblemBase::hasVariable(const std::string & var_name)
3405 {
3406  if (_nl->hasVariable(var_name))
3407  return true;
3408  else if (_aux->hasVariable(var_name))
3409  return true;
3410  else
3411  return false;
3412 }
3413 
3414 MooseVariable &
3415 FEProblemBase::getVariable(THREAD_ID tid, const std::string & var_name)
3416 {
3417  if (_nl->hasVariable(var_name))
3418  return _nl->getVariable(tid, var_name);
3419  else if (!_aux->hasVariable(var_name))
3420  mooseError("Unknown variable " + var_name);
3421 
3422  return _aux->getVariable(tid, var_name);
3423 }
3424 
3425 bool
3426 FEProblemBase::hasScalarVariable(const std::string & var_name)
3427 {
3428  if (_nl->hasScalarVariable(var_name))
3429  return true;
3430  else if (_aux->hasScalarVariable(var_name))
3431  return true;
3432  else
3433  return false;
3434 }
3435 
3437 FEProblemBase::getScalarVariable(THREAD_ID tid, const std::string & var_name)
3438 {
3439  if (_nl->hasScalarVariable(var_name))
3440  return _nl->getScalarVariable(tid, var_name);
3441  else if (_aux->hasScalarVariable(var_name))
3442  return _aux->getScalarVariable(tid, var_name);
3443  else
3444  mooseError("Unknown variable " + var_name);
3445 }
3446 
3447 System &
3448 FEProblemBase::getSystem(const std::string & var_name)
3449 {
3450  if (_nl->hasVariable(var_name))
3451  return _nl->system();
3452  else if (_aux->hasVariable(var_name))
3453  return _aux->system();
3454  else
3455  mooseError("Unable to find a system containing the variable " + var_name);
3456 }
3457 
3458 void
3459 FEProblemBase::setActiveElementalMooseVariables(const std::set<MooseVariable *> & moose_vars,
3460  THREAD_ID tid)
3461 {
3463 
3464  if (_displaced_problem)
3465  _displaced_problem->setActiveElementalMooseVariables(moose_vars, tid);
3466 }
3467 
3468 const std::set<MooseVariable *> &
3470 {
3472 }
3473 
3474 bool
3476 {
3478 }
3479 
3480 void
3482 {
3484 
3485  if (_displaced_problem)
3486  _displaced_problem->clearActiveElementalMooseVariables(tid);
3487 }
3488 
3489 void
3490 FEProblemBase::setActiveMaterialProperties(const std::set<unsigned int> & mat_prop_ids,
3491  THREAD_ID tid)
3492 {
3493  SubProblem::setActiveMaterialProperties(mat_prop_ids, tid);
3494 
3495  if (_displaced_problem)
3496  _displaced_problem->setActiveMaterialProperties(mat_prop_ids, tid);
3497 }
3498 
3499 const std::set<unsigned int> &
3501 {
3503 }
3504 
3505 bool
3507 {
3509 }
3510 
3511 void
3513 {
3515 
3516  if (_displaced_problem)
3517  _displaced_problem->clearActiveMaterialProperties(tid);
3518 }
3519 
3520 void
3521 FEProblemBase::createQRules(QuadratureType type, Order order, Order volume_order, Order face_order)
3522 {
3523  if (order == INVALID_ORDER)
3524  {
3525  // automatically determine the integration order
3526  order = _nl->getMinQuadratureOrder();
3527  if (order < _aux->getMinQuadratureOrder())
3528  order = _aux->getMinQuadratureOrder();
3529  }
3530 
3531  if (volume_order == INVALID_ORDER)
3532  volume_order = order;
3533 
3534  if (face_order == INVALID_ORDER)
3535  face_order = order;
3536 
3537  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
3538  _assembly[tid]->createQRules(type, order, volume_order, face_order);
3539 
3540  if (_displaced_problem)
3541  _displaced_problem->createQRules(type, order, volume_order, face_order);
3542 
3543  // Find the maximum number of quadrature points
3544  {
3545  MaxQpsThread mqt(*this, type, std::max(order, volume_order), face_order);
3546  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mqt);
3547  _max_qps = mqt.max();
3549 
3550  // If we have more shape functions or more quadrature points on
3551  // another processor, then we may need to handle those elements
3552  // ourselves later after repartitioning.
3553  _communicator.max(_max_qps);
3554  _communicator.max(_max_shape_funcs);
3555  }
3556 
3557  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
3558  {
3559  _zero[tid].resize(getMaxQps(), 0);
3560  _grad_zero[tid].resize(getMaxQps(), RealGradient(0.));
3561  _second_zero[tid].resize(getMaxQps(), RealTensor(0.));
3562  _second_phi_zero[tid].resize(getMaxQps(),
3563  std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));
3564  }
3565 }
3566 
3567 void
3569 {
3570  _coupling = type;
3571 }
3572 
3573 void
3575 {
3576  // TODO: Deprecate method
3577 
3579  _cm.reset(cm);
3580 }
3581 
3582 void
3583 FEProblemBase::setCouplingMatrix(std::unique_ptr<CouplingMatrix> cm)
3584 {
3586  _cm = std::move(cm);
3587 }
3588 
3589 void
3591 {
3592  unsigned int n_vars = _nl->nVariables();
3593  _nonlocal_cm.resize(n_vars);
3594  const auto & vars = _nl->getVariables(0);
3595  const auto & nonlocal_kernel = _nonlocal_kernels.getObjects();
3596  const auto & nonlocal_integrated_bc = _nonlocal_integrated_bcs.getObjects();
3597  for (const auto & ivar : vars)
3598  {
3599  for (const auto & kernel : nonlocal_kernel)
3600  {
3601  unsigned int i = ivar->number();
3602  if (i == kernel->variable().number())
3603  for (const auto & jvar : vars)
3604  {
3605  const auto it = _var_dof_map.find(jvar->name());
3606  if (it != _var_dof_map.end())
3607  {
3608  unsigned int j = jvar->number();
3609  _nonlocal_cm(i, j) = 1;
3610  }
3611  }
3612  }
3613  for (const auto & integrated_bc : nonlocal_integrated_bc)
3614  {
3615  unsigned int i = ivar->number();
3616  if (i == integrated_bc->variable().number())
3617  for (const auto & jvar : vars)
3618  {
3619  const auto it = _var_dof_map.find(jvar->name());
3620  if (it != _var_dof_map.end())
3621  {
3622  unsigned int j = jvar->number();
3623  _nonlocal_cm(i, j) = 1;
3624  }
3625  }
3626  }
3627  }
3628 }
3629 
3630 bool
3631 FEProblemBase::areCoupled(unsigned int ivar, unsigned int jvar)
3632 {
3633  return (*_cm)(ivar, jvar);
3634 }
3635 
3636 std::vector<std::pair<MooseVariable *, MooseVariable *>> &
3638 {
3639  return _assembly[tid]->couplingEntries();
3640 }
3641 
3642 std::vector<std::pair<MooseVariable *, MooseVariable *>> &
3644 {
3645  return _assembly[tid]->nonlocalCouplingEntries();
3646 }
3647 
3648 void
3650 {
3651  if (fe_cache)
3652  _console << "\nUtilizing FE Shape Function Caching\n" << std::endl;
3653 
3654  unsigned int n_threads = libMesh::n_threads();
3655 
3656  for (unsigned int i = 0; i < n_threads; ++i)
3657  _assembly[i]->useFECache(fe_cache); // fe_cache);
3658 }
3659 
3660 void
3662 {
3663  if (_initialized)
3664  return;
3665 
3666  unsigned int n_vars = _nl->nVariables();
3667  switch (_coupling)
3668  {
3669  case Moose::COUPLING_DIAG:
3670  _cm = libmesh_make_unique<CouplingMatrix>(n_vars);
3671  for (unsigned int i = 0; i < n_vars; i++)
3672  for (unsigned int j = 0; j < n_vars; j++)
3673  (*_cm)(i, j) = (i == j ? 1 : 0);
3674  break;
3675 
3676  // for full jacobian
3677  case Moose::COUPLING_FULL:
3678  _cm = libmesh_make_unique<CouplingMatrix>(n_vars);
3679  for (unsigned int i = 0; i < n_vars; i++)
3680  for (unsigned int j = 0; j < n_vars; j++)
3681  (*_cm)(i, j) = 1;
3682  break;
3683 
3685  // do nothing, _cm was already set through couplingMatrix() call
3686  break;
3687  }
3688 
3689  _nl->dofMap()._dof_coupling = _cm.get();
3690  _nl->dofMap().attach_extra_sparsity_function(&extraSparsity, _nl.get());
3691  _nl->dofMap().attach_extra_send_list_function(&extraSendList, _nl.get());
3692  _aux->dofMap().attach_extra_send_list_function(&extraSendList, _aux.get());
3693 
3694  if (_solve && n_vars == 0)
3695  mooseError("No variables specified in the FEProblemBase '", name(), "'.");
3696 
3697  ghostGhostedBoundaries(); // We do this again right here in case new boundaries have been added
3698 
3699  // do not assemble system matrix for JFNK solve
3700  if (solverParams()._type == Moose::ST_JFNK)
3701  _nl->turnOffJacobian();
3702 
3703  Moose::perf_log.push("eq.init()", "Setup");
3704  _eq.init();
3705  Moose::perf_log.pop("eq.init()", "Setup");
3706 
3707  Moose::perf_log.push("FEProblemBase::init::meshChanged()", "Setup");
3708  _mesh.meshChanged();
3709  if (_displaced_problem)
3711  Moose::perf_log.pop("FEProblemBase::init::meshChanged()", "Setup");
3712 
3713  Moose::perf_log.push("NonlinearSystem::update()", "Setup");
3714  _nl->update();
3715  Moose::perf_log.pop("NonlinearSystem::update()", "Setup");
3716 
3717  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3718  _assembly[tid]->init(_cm.get());
3719 
3720  _nl->init();
3721 
3722  if (_displaced_problem)
3723  _displaced_problem->init();
3724 
3725  _aux->init();
3726 
3727  _initialized = true;
3728 }
3729 
3730 void
3732 {
3733  Moose::perf_log.push("solve()", "Execution");
3734 
3735 #ifdef LIBMESH_HAVE_PETSC
3736  Moose::PetscSupport::petscSetOptions(*this); // Make sure the PETSc options are setup for this app
3737 #endif
3738 
3739  Moose::setSolverDefaults(*this);
3740 
3741  // Setup the output system for printing linear/nonlinear iteration information
3742  initPetscOutput();
3743 
3745 
3746  // reset flag so that linear solver does not use
3747  // the old converged reason "DIVERGED_NANORINF", when
3748  // we throw an exception and stop solve
3750 
3751  if (_solve)
3752  _nl->solve();
3753 
3754  if (_solve)
3755  _nl->update();
3756 
3757  // sync solutions in displaced problem
3758  if (_displaced_problem)
3759  _displaced_problem->syncSolutions();
3760 
3761  Moose::perf_log.pop("solve()", "Execution");
3762 }
3763 
3764 void
3765 FEProblemBase::setException(const std::string & message)
3766 {
3767  _has_exception = true;
3768  _exception_message = message;
3769 }
3770 
3771 void
3773 {
3774  // See if any processor had an exception. If it did, get back the
3775  // processor that the exception occurred on.
3776  unsigned int processor_id;
3777 
3778  _communicator.maxloc(_has_exception, processor_id);
3779 
3780  if (_has_exception)
3781  {
3782  _communicator.broadcast(_exception_message, processor_id);
3783 
3784  // Print the message
3785  if (_communicator.rank() == 0)
3786  Moose::err << _exception_message << std::endl;
3787 
3788  // Stop the solve -- this entails setting
3789  // SNESSetFunctionDomainError() or directly inserting NaNs in the
3790  // residual vector to let PETSc >= 3.6 return DIVERGED_NANORINF.
3791  _nl->stopSolve();
3792 
3793  // We've handled this exception, so we no longer have one.
3794  _has_exception = false;
3795 
3796  // Force the next linear convergence check to fail.
3798 
3799  // Repropagate the exception, so it can be caught at a higher level, typically
3800  // this is NonlinearSystem::computeResidual().
3802  }
3803 }
3804 
3805 bool
3807 {
3808  if (_solve)
3809  return _nl->converged();
3810  else
3811  return true;
3812 }
3813 
3814 unsigned int
3816 {
3817  return _nl->nNonlinearIterations();
3818 }
3819 
3820 unsigned int
3822 {
3823  return _nl->nLinearIterations();
3824 }
3825 
3826 Real
3828 {
3829  return _nl->finalNonlinearResidual();
3830 }
3831 
3832 bool
3834 {
3835  return _nl->computingInitialResidual();
3836 }
3837 
3838 void
3840 {
3841  _nl->copySolutionsBackwards();
3842  _aux->copySolutionsBackwards();
3843 }
3844 
3845 void
3847 {
3848  _nl->copyOldSolutions();
3849  _aux->copyOldSolutions();
3850 
3851  if (_displaced_problem != NULL)
3852  {
3853  _displaced_problem->nlSys().copyOldSolutions();
3854  _displaced_problem->auxSys().copyOldSolutions();
3855  }
3856 
3858 
3861 
3864 }
3865 
3866 void
3868 {
3869  _nl->restoreSolutions();
3870  _aux->restoreSolutions();
3871 
3872  if (_displaced_problem != NULL)
3873  _displaced_problem->updateMesh();
3874 }
3875 
3876 void
3878 {
3879  _nl->saveOldSolutions();
3880  _aux->saveOldSolutions();
3881 }
3882 
3883 void
3885 {
3886  _nl->restoreOldSolutions();
3887  _aux->restoreOldSolutions();
3888 }
3889 
3890 void
3892 {
3893  _nl->update();
3894  _aux->update();
3895  if (_displaced_problem != NULL)
3896  _displaced_problem->syncSolutions();
3898 }
3899 
3900 void
3902 {
3904 }
3905 
3906 void
3908 {
3910 }
3911 
3912 void
3914 {
3917 }
3918 
3919 void
3921 {
3922  _nl->onTimestepBegin();
3923 }
3924 
3925 void
3927 {
3928 }
3929 
3930 void
3932  const std::string & name,
3934 {
3935  setInputParametersFEProblem(parameters);
3936  parameters.set<SubProblem *>("_subproblem") = this;
3937  _aux->addTimeIntegrator(type, name + ":aux", parameters);
3938  _nl->addTimeIntegrator(type, name, parameters);
3939  _has_time_integrator = true;
3940 }
3941 
3942 void
3943 FEProblemBase::addPredictor(const std::string & type,
3944  const std::string & name,
3946 {
3947  setInputParametersFEProblem(parameters);
3948  parameters.set<SubProblem *>("_subproblem") = this;
3949  std::shared_ptr<Predictor> predictor = _factory.create<Predictor>(type, name, parameters);
3950  _nl->setPredictor(predictor);
3951 }
3952 
3953 Real
3955 {
3956  computeResidualType(*_nl->currentSolution(), _nl->RHS());
3957 
3958  return _nl->RHS().l2_norm();
3959 }
3960 
3961 void
3962 FEProblemBase::computeResidual(NonlinearImplicitSystem & /*sys*/,
3963  const NumericVector<Number> & soln,
3964  NumericVector<Number> & residual)
3965 {
3966  computeResidual(soln, residual);
3967 }
3968 
3969 void
3970 FEProblemBase::computeResidual(const NumericVector<Number> & soln, NumericVector<Number> & residual)
3971 {
3972  try
3973  {
3974  computeResidualType(soln, residual, _kernel_type);
3975  }
3976  catch (MooseException & e)
3977  {
3978  // If a MooseException propagates all the way to here, it means
3979  // that it was thrown from a MOOSE system where we do not
3980  // (currently) properly support the throwing of exceptions, and
3981  // therefore we have no choice but to error out. It may be
3982  // *possible* to handle exceptions from other systems, but in the
3983  // meantime, we don't want to silently swallow any unhandled
3984  // exceptions here.
3985  mooseError("An unhandled MooseException was raised during residual computation. Please "
3986  "contact the MOOSE team for assistance.");
3987  }
3988 }
3989 
3990 void
3992  const NumericVector<Number> & u,
3993  const NumericVector<Number> & udot,
3994  NumericVector<Number> & residual)
3995 {
3996  _nl->setSolutionUDot(udot);
3997  _time = time;
3998  computeResidual(u, residual);
3999 }
4000 
4001 void
4002 FEProblemBase::computeResidualType(const NumericVector<Number> & soln,
4003  NumericVector<Number> & residual,
4005 {
4006  _nl->setSolution(soln);
4007 
4008  _nl->zeroVariablesForResidual();
4009  _aux->zeroVariablesForResidual();
4010 
4011  unsigned int n_threads = libMesh::n_threads();
4012 
4014 
4015  // Random interface objects
4016  for (const auto & it : _random_data_objects)
4017  it.second->updateSeeds(EXEC_LINEAR);
4018 
4020 
4022 
4023  for (unsigned int tid = 0; tid < n_threads; tid++)
4024  reinitScalars(tid);
4025 
4027 
4028  if (_displaced_problem != NULL)
4029  _displaced_problem->updateMesh();
4030 
4031  for (THREAD_ID tid = 0; tid < n_threads; tid++)
4032  {
4035  }
4036  _aux->residualSetup();
4037 
4038  _nl->computeTimeDerivatives();
4039 
4040  try
4041  {
4042  _aux->compute(EXEC_LINEAR);
4043  }
4044  catch (MooseException & e)
4045  {
4046  _console << "\nA MooseException was raised during Auxiliary variable computation.\n"
4047  << "The next solve will fail, the timestep will be reduced, and we will try again.\n"
4048  << std::endl;
4049 
4050  // We know the next solve is going to fail, so there's no point in
4051  // computing anything else after this. Plus, using incompletely
4052  // computed AuxVariables in subsequent calculations could lead to
4053  // other errors or unhandled exceptions being thrown.
4054  return;
4055  }
4056 
4058 
4060 
4062 
4064 
4065  _nl->computeResidual(residual, type);
4066 }
4067 
4068 void
4069 FEProblemBase::computeJacobian(NonlinearImplicitSystem & /*sys*/,
4070  const NumericVector<Number> & soln,
4071  SparseMatrix<Number> & jacobian)
4072 {
4073  computeJacobian(soln, jacobian, Moose::KT_ALL);
4074 }
4075 
4076 void
4077 FEProblemBase::computeJacobian(const NumericVector<Number> & soln,
4078  SparseMatrix<Number> & jacobian,
4079  Moose::KernelType kernel_type)
4080 {
4081  if (!_has_jacobian || !_const_jacobian)
4082  {
4083  _nl->setSolution(soln);
4084 
4085  _nl->zeroVariablesForJacobian();
4086  _aux->zeroVariablesForJacobian();
4087 
4088  unsigned int n_threads = libMesh::n_threads();
4089 
4090  // Random interface objects
4091  for (const auto & it : _random_data_objects)
4092  it.second->updateSeeds(EXEC_NONLINEAR);
4093 
4096 
4099 
4100  for (unsigned int tid = 0; tid < n_threads; tid++)
4101  reinitScalars(tid);
4102 
4104 
4105  if (_displaced_problem != NULL)
4106  _displaced_problem->updateMesh();
4107 
4108  for (unsigned int tid = 0; tid < n_threads; tid++)
4109  {
4112  }
4113 
4114  _aux->jacobianSetup();
4115 
4116  _aux->compute(EXEC_NONLINEAR);
4117 
4119 
4121 
4123 
4124  _nl->computeJacobian(jacobian, kernel_type);
4125 
4128  _has_jacobian = true;
4129  }
4130 
4132  {
4133  // This call is here to make sure the residual vector is up to date with any decisions that have
4134  // been made in
4135  // the Jacobian evaluation. That is important in JFNK because that residual is used for finite
4136  // differencing
4137  computeResidual(soln, _nl->RHS());
4138  _nl->RHS().close();
4139  }
4140 }
4141 
4142 void
4144  const NumericVector<Number> & u,
4145  const NumericVector<Number> & udot,
4146  Real shift,
4147  SparseMatrix<Number> & jacobian)
4148 {
4149  if (0)
4150  { // The current interface guarantees that the residual is called before Jacobian, thus udot has
4151  // already been set
4152  _nl->setSolutionUDot(udot);
4153  }
4154  _nl->duDotDu() = shift;
4155  _time = time;
4156  computeJacobian(u, jacobian);
4157 }
4158 
4159 void
4160 FEProblemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
4161 {
4162  if (_displaced_problem != NULL)
4163  _displaced_problem->updateMesh();
4164 
4165  _aux->compute(EXEC_NONLINEAR);
4166 
4167  _nl->computeJacobianBlocks(blocks);
4168 }
4169 
4170 void
4171 FEProblemBase::computeJacobianBlock(SparseMatrix<Number> & jacobian,
4172  libMesh::System & precond_system,
4173  unsigned int ivar,
4174  unsigned int jvar)
4175 {
4176  std::vector<JacobianBlock *> blocks;
4177  JacobianBlock * block = new JacobianBlock(precond_system, jacobian, ivar, jvar);
4178  blocks.push_back(block);
4179  computeJacobianBlocks(blocks);
4180  delete block;
4181 }
4182 
4183 void
4184 FEProblemBase::computeBounds(NonlinearImplicitSystem & /*sys*/,
4185  NumericVector<Number> & lower,
4186  NumericVector<Number> & upper)
4187 {
4188  if (!_nl->hasVector("lower_bound") || !_nl->hasVector("upper_bound"))
4189  return;
4190 
4191  NumericVector<Number> & _lower = _nl->getVector("lower_bound");
4192  NumericVector<Number> & _upper = _nl->getVector("upper_bound");
4193  _lower.swap(lower);
4194  _upper.swap(upper);
4195  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
4197 
4198  _aux->residualSetup();
4199  _aux->compute(EXEC_LINEAR);
4200  _lower.swap(lower);
4201  _upper.swap(upper);
4202 
4204 }
4205 
4206 void
4207 FEProblemBase::computeNearNullSpace(NonlinearImplicitSystem & /*sys*/,
4208  std::vector<NumericVector<Number> *> & sp)
4209 {
4210  sp.clear();
4211  for (unsigned int i = 0; i < subspaceDim("NearNullSpace"); ++i)
4212  {
4213  std::stringstream postfix;
4214  postfix << "_" << i;
4215  std::string modename = "NearNullSpace" + postfix.str();
4216  sp.push_back(&_nl->getVector(modename));
4217  }
4218 }
4219 
4220 void
4221 FEProblemBase::computeNullSpace(NonlinearImplicitSystem & /*sys*/,
4222  std::vector<NumericVector<Number> *> & sp)
4223 {
4224  sp.clear();
4225  for (unsigned int i = 0; i < subspaceDim("NullSpace"); ++i)
4226  {
4227  std::stringstream postfix;
4228  postfix << "_" << i;
4229  sp.push_back(&_nl->getVector("NullSpace" + postfix.str()));
4230  }
4231 }
4232 
4233 void
4234 FEProblemBase::computeTransposeNullSpace(NonlinearImplicitSystem & /*sys*/,
4235  std::vector<NumericVector<Number> *> & sp)
4236 {
4237  sp.clear();
4238  for (unsigned int i = 0; i < subspaceDim("TransposeNullSpace"); ++i)
4239  {
4240  std::stringstream postfix;
4241  postfix << "_" << i;
4242  sp.push_back(&_nl->getVector("TransposeNullSpace" + postfix.str()));
4243  }
4244 }
4245 
4246 void
4247 FEProblemBase::computePostCheck(NonlinearImplicitSystem & sys,
4248  const NumericVector<Number> & old_soln,
4249  NumericVector<Number> & search_direction,
4250  NumericVector<Number> & new_soln,
4251  bool & changed_search_direction,
4252  bool & changed_new_soln)
4253 {
4254  // This function replaces the old PetscSupport::dampedCheck() function.
4255  //
4256  // 1.) Recreate code in PetscSupport::dampedCheck() for constructing
4257  // ghosted "soln" and "update" vectors.
4258  // 2.) Call FEProblemBase::computeDamping() with these ghost vectors.
4259  // 3.) Recreate the code in PetscSupport::dampedCheck() to actually update
4260  // the solution vector based on the damping, and set the "changed" flags
4261  // appropriately.
4262  Moose::perf_log.push("computePostCheck()", "Execution");
4263 
4264  // MOOSE's FEProblemBase doesn't update the solution during the
4265  // postcheck, but FEProblemBase-derived classes might.
4267  {
4268  // We need ghosted versions of new_soln and search_direction (the
4269  // ones we get from libmesh/PETSc are PARALLEL vectors. To make
4270  // our lives simpler, we use the same ghosting pattern as the
4271  // system's current_local_solution to create new ghosted vectors.
4272 
4273  // Construct zeroed-out clones with the same ghosted dofs as the
4274  // System's current_local_solution.
4275  std::unique_ptr<NumericVector<Number>> ghosted_solution =
4276  sys.current_local_solution->zero_clone(),
4277  ghosted_search_direction =
4278  sys.current_local_solution->zero_clone();
4279 
4280  // Copy values from input vectors into clones with ghosted values.
4281  *ghosted_solution = new_soln;
4282  *ghosted_search_direction = search_direction;
4283 
4284  if (_has_dampers)
4285  {
4286  // Compute the damping coefficient using the ghosted vectors
4287  Real damping = computeDamping(*ghosted_solution, *ghosted_search_direction);
4288 
4289  // If some non-trivial damping was computed, update the new_soln
4290  // vector accordingly.
4291  if (damping < 1.0)
4292  {
4293  new_soln = old_soln;
4294  new_soln.add(-damping, search_direction);
4295  changed_new_soln = true;
4296  }
4297  }
4298 
4299  if (shouldUpdateSolution())
4300  {
4301  // Update the ghosted copy of the new solution, if necessary.
4302  if (changed_new_soln)
4303  *ghosted_solution = new_soln;
4304 
4305  bool updated_solution = updateSolution(new_soln, *ghosted_solution);
4306  if (updated_solution)
4307  changed_new_soln = true;
4308  }
4309  }
4310 
4312  {
4313  _nl->setPreviousNewtonSolution(old_soln);
4314  _aux->setPreviousNewtonSolution();
4315  }
4316 
4317  // MOOSE doesn't change the search_direction
4318  changed_search_direction = false;
4319 
4320  Moose::perf_log.pop("computePostCheck()", "Execution");
4321 }
4322 
4323 Real
4324 FEProblemBase::computeDamping(const NumericVector<Number> & soln,
4325  const NumericVector<Number> & update)
4326 {
4327  Moose::perf_log.push("compute_dampers()", "Execution");
4328 
4329  // Default to no damping
4330  Real damping = 1.0;
4331 
4332  if (_has_dampers)
4333  {
4334  // Save pointer to the current solution
4335  const NumericVector<Number> * _saved_current_solution = _nl->currentSolution();
4336 
4337  _nl->setSolution(soln);
4338  // For now, do not re-compute auxiliary variables. Doing so allows a wild solution increment
4339  // to get to the material models, which may not be able to cope with drastically different
4340  // values. Once more complete dependency checking is in place, auxiliary variables (and
4341  // material properties) will be computed as needed by dampers.
4342  // _aux.compute();
4343  damping = _nl->computeDamping(soln, update);
4344 
4345  // restore saved solution
4346  _nl->setSolution(*_saved_current_solution);
4347  }
4348 
4349  Moose::perf_log.pop("compute_dampers()", "Execution");
4350 
4351  return damping;
4352 }
4353 
4354 bool
4356 {
4357  return false;
4358 }
4359 
4360 bool
4361 FEProblemBase::updateSolution(NumericVector<Number> & /*vec_solution*/,
4362  NumericVector<Number> & /*ghosted_solution*/)
4363 {
4364  return false;
4365 }
4366 
4367 void
4368 FEProblemBase::predictorCleanup(NumericVector<Number> & /*ghosted_solution*/)
4369 {
4370 }
4371 
4372 void
4373 FEProblemBase::addDisplacedProblem(std::shared_ptr<DisplacedProblem> displaced_problem)
4374 {
4375  _displaced_mesh = &displaced_problem->mesh();
4376  _displaced_problem = displaced_problem;
4377 }
4378 
4379 void
4381 {
4383 
4384  if (_displaced_problem)
4385  _displaced_problem->updateGeomSearch(type);
4386 }
4387 
4388 void
4390 {
4391  if (_displaced_problem) // Only need to do this if things are moving...
4392  {
4393  switch (_mesh.getPatchUpdateStrategy())
4394  {
4395  case 0: // Never
4396  break;
4397  case 2: // Auto
4398  {
4399  Real max = _displaced_problem->geomSearchData().maxPatchPercentage();
4400  _communicator.max(max);
4401 
4402  // If we haven't moved very far through the patch
4403  if (max < 0.4)
4404  break;
4405  }
4406 
4407  // Let this fall through if things do need to be updated...
4408 
4409  case 1: // Always
4410  // Flush output here to see the message before the reinitialization, which could take a
4411  // while
4412  _console << "\n\nUpdating geometric search patches\n" << std::endl;
4413 
4416 
4417  _displaced_problem->geomSearchData().clearNearestNodeLocators();
4419 
4421 
4422  // This is needed to reinitialize PETSc output
4423  initPetscOutput();
4424  }
4425  }
4426 }
4427 
4428 #ifdef LIBMESH_ENABLE_AMR
4429 void
4431 {
4432  unsigned int n = adaptivity().getInitialSteps();
4433  _cycles_completed = 0;
4434  for (unsigned int i = 0; i < n; i++)
4435  {
4436  _console << "Initial adaptivity step " << i + 1 << " of " << n << std::endl;
4438  computeMarkers();
4439 
4441  {
4442  meshChanged();
4443 
4444  // reproject the initial condition
4445  projectSolution();
4446 
4448  }
4449  else
4450  {
4451  _console << "Mesh unchanged, skipping remaining steps..." << std::endl;
4452  return;
4453  }
4454  }
4455 }
4456 
4457 void
4459 {
4460  // reset cycle counter
4461  _cycles_completed = 0;
4462 
4464  return;
4465 
4466  unsigned int cycles_per_step = _adaptivity.getCyclesPerStep();
4467 
4468  Moose::perf_log.push("Adaptivity: adaptMesh()", "Execution");
4469 
4470  for (unsigned int i = 0; i < cycles_per_step; ++i)
4471  {
4472  _console << "Adaptivity step " << i + 1 << " of " << cycles_per_step << '\n';
4473  // Markers were already computed once by Executioner
4474  if (_adaptivity.getRecomputeMarkersFlag() && i > 0)
4475  computeMarkers();
4476  if (_adaptivity.adaptMesh())
4477  {
4478  meshChanged();
4480  }
4481  else
4482  {
4483  _console << "Mesh unchanged, skipping remaining steps..." << std::endl;
4484  break;
4485  }
4486 
4487  // Show adaptivity progress
4488  _console << std::flush;
4489  }
4490 
4491  Moose::perf_log.pop("Adaptivity: adaptMesh()", "Execution");
4492 }
4493 #endif // LIBMESH_ENABLE_AMR
4494 
4495 void
4496 FEProblemBase::initXFEM(std::shared_ptr<XFEMInterface> xfem)
4497 {
4498  _xfem = xfem;
4499  _xfem->setMesh(&_mesh);
4500  if (_displaced_mesh)
4501  _xfem->setDisplacedMesh(_displaced_mesh);
4502  _xfem->setMaterialData(&_material_data);
4503  _xfem->setBoundaryMaterialData(&_bnd_material_data);
4504 
4505  unsigned int n_threads = libMesh::n_threads();
4506  for (unsigned int i = 0; i < n_threads; ++i)
4507  {
4508  _assembly[i]->setXFEM(_xfem);
4509  if (_displaced_problem != NULL)
4510  _displaced_problem->assembly(i).setXFEM(_xfem);
4511  }
4512 }
4513 
4514 bool
4516 {
4517  bool updated = false;
4518  if (haveXFEM())
4519  {
4520  updated = _xfem->update(_time, *_nl, *_aux);
4521  if (updated)
4522  {
4523  meshChanged();
4524  _xfem->initSolution(*_nl, *_aux);
4525  restoreSolutions();
4526  }
4527  }
4528  return updated;
4529 }
4530 
4531 void
4533 {
4535  _mesh.cacheChangedLists(); // Currently only used with adaptivity and stateful material
4536  // properties
4537 
4538  // Clear these out because they corresponded to the old mesh
4539  _ghosted_elems.clear();
4540 
4542 
4543  // The mesh changed. We notify the MooseMesh first, because
4544  // callbacks (e.g. for sparsity calculations) triggered by the
4545  // EquationSystems reinit may require up-to-date MooseMesh caches.
4546  _mesh.meshChanged();
4547  _eq.reinit();
4548 
4549  // But that breaks other adaptivity code, unless we then *again*
4550  // update the MooseMesh caches. E.g. the definition of "active" and
4551  // "local" may be *changed* by EquationSystems::reinit().
4552  _mesh.meshChanged();
4553 
4554  // Since the Mesh changed, update the PointLocator object used by DiracKernels.
4556 
4557  unsigned int n_threads = libMesh::n_threads();
4558 
4559  for (unsigned int i = 0; i < n_threads; ++i)
4560  _assembly[i]->invalidateCache();
4561 
4562  // Need to redo ghosting
4564 
4565  if (_displaced_problem != NULL)
4566  {
4567  _displaced_problem->meshChanged();
4569  }
4570 
4572 
4574 
4575  // We need to create new storage for the new elements and copy stateful properties from the old
4576  // elements.
4579  {
4580  {
4581  ProjectMaterialProperties pmp(true,
4582  *this,
4583  *_nl,
4588  _assembly);
4589  Threads::parallel_reduce(*_mesh.refinedElementRange(), pmp);
4590  }
4591 
4592  {
4593  ProjectMaterialProperties pmp(false,
4594  *this,
4595  *_nl,
4600  _assembly);
4601  Threads::parallel_reduce(*_mesh.coarsenedElementRange(), pmp);
4602  }
4603  }
4604 
4607 
4608  _has_jacobian = false; // we have to recompute jacobian when mesh changed
4609 
4610  for (const auto & mci : _notify_when_mesh_changes)
4611  mci->meshChanged();
4612 }
4613 
4614 void
4616 {
4617  _notify_when_mesh_changes.push_back(mci);
4618 }
4619 
4620 void
4622 {
4623  // Check for unsatisfied actions
4624  const std::set<SubdomainID> & mesh_subdomains = _mesh.meshSubdomains();
4625 
4626  // Check kernel coverage of subdomains (blocks) in the mesh
4628  _nl->checkKernelCoverage(mesh_subdomains);
4629 
4630  // Check materials
4631  {
4632 #ifdef LIBMESH_ENABLE_AMR
4633  if (_adaptivity.isOn() &&
4635  {
4636  _console << "Using EXPERIMENTAL Stateful Material Property projection with Adaptivity!\n";
4637 
4638  if (n_processors() > 1)
4639  {
4640  if (_mesh.uniformRefineLevel() > 0 && _mesh.getMesh().skip_partitioning() == false)
4641  mooseError("This simulation is using uniform refinement on the mesh, with stateful "
4642  "properties and adaptivity. "
4643  "You must skip partitioning to run this case:\nMesh/skip_partitioning=true");
4644 
4645  _console << "\nWarning! Mesh re-partitioning is disabled while using stateful material "
4646  "properties! This can lead to large load imbalances and degraded "
4647  "performance!!\n\n";
4648  _mesh.getMesh().skip_partitioning(true);
4649 
4650  _mesh.errorIfDistributedMesh("StatefulMaterials + Adaptivity");
4651 
4652  if (_displaced_problem)
4653  _displaced_problem->mesh().getMesh().skip_partitioning(true);
4654  }
4655  }
4656 #endif
4657 
4658  std::set<SubdomainID> local_mesh_subs(mesh_subdomains);
4659 
4661  {
4666  bool check_material_coverage = false;
4667  std::set<SubdomainID> ids = _all_materials.getActiveBlocks();
4668  for (const auto & id : ids)
4669  {
4670  local_mesh_subs.erase(id);
4671  check_material_coverage = true;
4672  }
4673 
4674  // also exclude mortar spaces from the material check
4675  auto & mortar_ifaces = _mesh.getMortarInterfaces();
4676  for (const auto & mortar_iface : mortar_ifaces)
4677  local_mesh_subs.erase(mortar_iface->_id);
4678 
4679  // Check Material Coverage
4680  if (check_material_coverage && !local_mesh_subs.empty())
4681  {
4682  std::stringstream extra_subdomain_ids;
4684  std::copy(local_mesh_subs.begin(),
4685  local_mesh_subs.end(),
4686  std::ostream_iterator<unsigned int>(extra_subdomain_ids, " "));
4687 
4688  mooseError("The following blocks from your input mesh do not contain an active material: " +
4689  extra_subdomain_ids.str() + "\nWhen ANY mesh block contains a Material object, "
4690  "all blocks must contain a Material object.\n");
4691  }
4692  }
4693 
4694  // Check material properties on blocks and boundaries
4697 
4698  // Check that material properties exist when requested by other properties on a given block
4699  const auto & materials = _all_materials.getActiveObjects();
4700  for (const auto & material : materials)
4701  material->checkStatefulSanity();
4702 
4704  }
4705 
4706  // Check UserObjects and Postprocessors
4707  checkUserObjects();
4708 
4709  // Verify that we don't have any Element type/Coordinate Type conflicts
4711 
4712  // If using displacements, verify that the order of the displacement
4713  // variables matches the order of the elements in the displaced
4714  // mesh.
4716 }
4717 
4718 void
4720 {
4721  if (_displaced_problem)
4722  {
4723  MeshBase::const_element_iterator it = _displaced_mesh->activeLocalElementsBegin(),
4725 
4726  bool mesh_has_second_order_elements = false;
4727  for (; it != end; ++it)
4728  {
4729  if ((*it)->default_order() == SECOND)
4730  {
4731  mesh_has_second_order_elements = true;
4732  break;
4733  }
4734  }
4735 
4736  // We checked our local elements, so take the max over all processors.
4737  _displaced_mesh->comm().max(mesh_has_second_order_elements);
4738 
4739  // If the Mesh has second order elements, make sure the
4740  // displacement variables are second-order.
4741  if (mesh_has_second_order_elements)
4742  {
4743  const std::vector<std::string> & displacement_variables =
4744  _displaced_problem->getDisplacementVarNames();
4745 
4746  for (const auto & var_name : displacement_variables)
4747  {
4748  MooseVariable & mv = _displaced_problem->getVariable(/*tid=*/0, var_name);
4749  if (mv.order() != SECOND)
4750  mooseError("Error: mesh has SECOND order elements, so all displacement variables must be "
4751  "SECOND order.");
4752  }
4753  }
4754  }
4755 }
4756 
4757 void
4759 {
4760  // Check user_objects block coverage
4761  std::set<SubdomainID> mesh_subdomains = _mesh.meshSubdomains();
4762  std::set<SubdomainID> user_objects_blocks;
4763 
4764  // gather names of all user_objects that were defined in the input file
4765  // and the blocks that they are defined on
4766  std::set<std::string> names;
4767 
4768  const auto & objects = _all_user_objects.getActiveObjects();
4769  for (const auto & obj : objects)
4770  names.insert(obj->name());
4771 
4772  // See if all referenced blocks are covered
4773  mesh_subdomains.insert(Moose::ANY_BLOCK_ID);
4774  std::set<SubdomainID> difference;
4775  std::set_difference(user_objects_blocks.begin(),
4776  user_objects_blocks.end(),
4777  mesh_subdomains.begin(),
4778  mesh_subdomains.end(),
4779  std::inserter(difference, difference.end()));
4780 
4781  if (!difference.empty())
4782  {
4783  std::ostringstream oss;
4784  oss << "One or more UserObjects is referencing a nonexistent block:\n";
4785  for (const auto & id : difference)
4786  oss << id << "\n";
4787  mooseError(oss.str());
4788  }
4789 
4790  // check that all requested UserObjects were defined in the input file
4791  for (const auto & it : _pps_data.values())
4792  {
4793  if (names.find(it.first) == names.end())
4794  mooseError("Postprocessor '" + it.first + "' requested but not specified in the input file.");
4795  }
4796 }
4797 
4798 void
4800  const std::map<SubdomainID, std::vector<std::shared_ptr<Material>>> & materials_map)
4801 {
4802  auto & prop_names = _material_props.statefulPropNames();
4803 
4804  for (const auto & it : materials_map)
4805  {
4807  std::set<std::string> block_depend_props, block_supplied_props;
4808 
4809  for (const auto & mat1 : it.second)
4810  {
4811  const std::set<std::string> & depend_props = mat1->getRequestedItems();
4812  block_depend_props.insert(depend_props.begin(), depend_props.end());
4813 
4814  auto & alldeps = mat1->getMatPropDependencies(); // includes requested stateful props
4815  for (auto & dep : alldeps)
4816  {
4817  if (prop_names.count(dep) > 0)
4818  block_depend_props.insert(prop_names.at(dep));
4819  }
4820 
4821  // See if any of the active materials supply this property
4822  for (const auto & mat2 : it.second)
4823  {
4824  const std::set<std::string> & supplied_props = mat2->Material::getSuppliedItems();
4825  block_supplied_props.insert(supplied_props.begin(), supplied_props.end());
4826  }
4827  }
4828 
4829  // Add zero material properties specific to this block and unrestricted
4830  block_supplied_props.insert(_zero_block_material_props[it.first].begin(),
4831  _zero_block_material_props[it.first].end());
4832  block_supplied_props.insert(_zero_block_material_props[Moose::ANY_BLOCK_ID].begin(),
4834 
4835  // Error check to make sure all properties consumed by materials are supplied on this block
4836  std::set<std::string> difference;
4837  std::set_difference(block_depend_props.begin(),
4838  block_depend_props.end(),
4839  block_supplied_props.begin(),
4840  block_supplied_props.end(),
4841  std::inserter(difference, difference.end()));
4842 
4843  if (!difference.empty())
4844  {
4845  std::ostringstream oss;
4846  oss << "One or more Material Properties were not supplied on block " << it.first << ":\n";
4847  for (const auto & name : difference)
4848  oss << name << "\n";
4849  mooseError(oss.str());
4850  }
4851  }
4852 
4853  // This loop checks that materials are not supplied by multiple Material objects
4854  for (const auto & it : materials_map)
4855  {
4856  const auto & materials = it.second;
4857  std::set<std::string> inner_supplied, outer_supplied;
4858 
4859  for (const auto & outer_mat : materials)
4860  {
4861  // Storage for properties for this material (outer) and all other materials (inner)
4862  outer_supplied = outer_mat->getSuppliedItems();
4863  inner_supplied.clear();
4864 
4865  // Property to material map for error reporting
4866  std::map<std::string, std::set<std::string>> prop_to_mat;
4867  for (const auto & name : outer_supplied)
4868  prop_to_mat[name].insert(outer_mat->name());
4869 
4870  for (const auto & inner_mat : materials)
4871  {
4872  if (outer_mat == inner_mat)
4873  continue;
4874  inner_supplied.insert(inner_mat->getSuppliedItems().begin(),
4875  inner_mat->getSuppliedItems().end());
4876 
4877  for (const auto & inner_supplied_name : inner_supplied)
4878  prop_to_mat[inner_supplied_name].insert(inner_mat->name());
4879  }
4880 
4881  // Test that a property isn't supplied on multiple blocks
4882  std::set<std::string> intersection;
4883  std::set_intersection(outer_supplied.begin(),
4884  outer_supplied.end(),
4885  inner_supplied.begin(),
4886  inner_supplied.end(),
4887  std::inserter(intersection, intersection.end()));
4888 
4889  if (!intersection.empty())
4890  {
4891  std::ostringstream oss;
4892  oss << "The following material properties are declared on block " << it.first
4893  << " by multiple materials:\n";
4894  oss << ConsoleUtils::indent(2) << std::setw(30) << std::left << "Material Property"
4895  << "Material Objects\n";
4896  for (const auto & outer_name : intersection)
4897  {
4898  oss << ConsoleUtils::indent(2) << std::setw(30) << std::left << outer_name;
4899  for (const auto & inner_name : prop_to_mat[outer_name])
4900  oss << inner_name << " ";
4901  oss << '\n';
4902  }
4903 
4904  mooseError(oss.str());
4905  break;
4906  }
4907  }
4908  }
4909 }
4910 
4911 void
4913 {
4914  for (const auto & elem : _mesh.getMesh().element_ptr_range())
4915  {
4916  SubdomainID sid = elem->subdomain_id();
4917  if (_coord_sys[sid] == Moose::COORD_RZ && elem->dim() == 3)
4918  mooseError("An RZ coordinate system was requested for subdomain " + Moose::stringify(sid) +
4919  " which contains 3D elements.");
4920  if (_coord_sys[sid] == Moose::COORD_RSPHERICAL && elem->dim() > 1)
4921  mooseError("An RSPHERICAL coordinate system was requested for subdomain " +
4922  Moose::stringify(sid) + " which contains 2D or 3D elements.");
4923  }
4924 }
4925 
4926 void
4928 {
4929  _nl->serializeSolution();
4930  _aux->serializeSolution();
4931 }
4932 
4933 void
4934 FEProblemBase::setRestartFile(const std::string & file_name)
4935 {
4936  _app.setRestart(true);
4937  _resurrector->setRestartFile(file_name);
4938  if (_app.getRecoverFileSuffix() == "cpa")
4939  _resurrector->setRestartSuffix("xda");
4940 }
4941 
4942 std::vector<VariableName>
4944 {
4945  std::vector<VariableName> names;
4946 
4947  const std::vector<VariableName> & nl_var_names = _nl->getVariableNames();
4948  names.insert(names.end(), nl_var_names.begin(), nl_var_names.end());
4949 
4950  const std::vector<VariableName> & aux_var_names = _aux->getVariableNames();
4951  names.insert(names.end(), aux_var_names.begin(), aux_var_names.end());
4952 
4953  return names;
4954 }
4955 
4958  const PetscInt it,
4959  const Real xnorm,
4960  const Real snorm,
4961  const Real fnorm,
4962  const Real rtol,
4963  const Real stol,
4964  const Real abstol,
4965  const PetscInt nfuncs,
4966  const PetscInt max_funcs,
4967  const Real initial_residual_before_preset_bcs,
4968  const Real div_threshold)
4969 {
4972 
4973  // This is the first residual before any iterations have been done,
4974  // but after PresetBCs (if any) have been imposed on the solution
4975  // vector. We save it, and use it to detect convergence if
4976  // compute_initial_residual_before_preset_bcs=false.
4977  if (it == 0)
4978  system._initial_residual_after_preset_bcs = fnorm;
4979 
4980  std::ostringstream oss;
4981  if (fnorm != fnorm)
4982  {
4983  oss << "Failed to converge, function norm is NaN\n";
4984  reason = MOOSE_DIVERGED_FNORM_NAN;
4985  }
4986  else if (fnorm < abstol)
4987  {
4988  oss << "Converged due to function norm " << fnorm << " < " << abstol << '\n';
4989  reason = MOOSE_CONVERGED_FNORM_ABS;
4990  }
4991  else if (nfuncs >= max_funcs)
4992  {
4993  oss << "Exceeded maximum number of function evaluations: " << nfuncs << " > " << max_funcs
4994  << '\n';
4996  }
4997  else if (it && fnorm > system._last_nl_rnorm && fnorm >= div_threshold)
4998  {
4999  oss << "Nonlinear solve was blowing up!\n";
5000  reason = MOOSE_DIVERGED_LINE_SEARCH;
5001  }
5002 
5003  if (it && !reason)
5004  {
5005  // If compute_initial_residual_before_preset_bcs==false, then use the
5006  // first residual computed by Petsc to determine convergence.
5007  Real the_residual = system._compute_initial_residual_before_preset_bcs
5008  ? initial_residual_before_preset_bcs
5010  if (fnorm <= the_residual * rtol)
5011  {
5012  oss << "Converged due to function norm " << fnorm << " < "
5013  << " (relative tolerance)\n";
5015  }
5016  else if (snorm < stol * xnorm)
5017  {
5018  oss << "Converged due to small update length: " << snorm << " < " << stol << " * " << xnorm
5019  << '\n';
5021  }
5022  }
5023 
5024  system._last_nl_rnorm = fnorm;
5025  system._current_nl_its = static_cast<unsigned int>(it);
5026 
5027  msg = oss.str();
5028  if (_app.multiAppLevel() > 0)
5030 
5031  return (reason);
5032 }
5033 
5036  const PetscInt n,
5037  const Real rnorm,
5038  const Real /*rtol*/,
5039  const Real /*atol*/,
5040  const Real /*dtol*/,
5041  const PetscInt maxits)
5042 {
5044  {
5045  // Unset the flag
5047  return MOOSE_DIVERGED_NANORINF;
5048  }
5049 
5050  // We initialize the reason to something that basically means MOOSE
5051  // has not made a decision on convergence yet.
5053 
5054  // Get a reference to our Nonlinear System
5056 
5057  // If it's the beginning of a new set of iterations, reset
5058  // last_rnorm, otherwise record the most recent linear residual norm
5059  // in the NonlinearSystem.
5060  if (n == 0)
5061  system._last_rnorm = 1e99;
5062  else
5063  system._last_rnorm = rnorm;
5064 
5065  // If the linear residual norm is less than the System's linear absolute
5066  // step tolerance, we consider it to be converged and set the reason as
5067  // MOOSE_CONVERGED_RTOL.
5068  if (std::abs(rnorm - system._last_rnorm) < system._l_abs_step_tol)
5069  reason = MOOSE_CONVERGED_RTOL;
5070 
5071  // If we hit max its, then we consider that converged (rather than
5072  // KSP_DIVERGED_ITS).
5073  if (n >= maxits)
5074  reason = MOOSE_CONVERGED_ITS;
5075 
5076  // If either of our convergence criteria is met, store the number of linear
5077  // iterations in the System.
5078  if (reason == MOOSE_CONVERGED_ITS || reason == MOOSE_CONVERGED_RTOL)
5079  system._current_l_its.push_back(static_cast<unsigned int>(n));
5080 
5081  return reason;
5082 }
5083 
5084 SolverParams &
5086 {
5087  return _solver_params;
5088 }
5089 
5090 void
5091 FEProblemBase::registerRandomInterface(RandomInterface & random_interface, const std::string & name)
5092 {
5093  auto insert_pair = moose_try_emplace(
5094  _random_data_objects, name, libmesh_make_unique<RandomData>(*this, random_interface));
5095 
5096  auto random_data_ptr = insert_pair.first->second.get();
5097  random_interface.setRandomDataPointer(random_data_ptr);
5098 }
5099 
5100 bool
5102 {
5103  if (_bnd_mat_side_cache[tid].find(bnd_id) == _bnd_mat_side_cache[tid].end())
5104  {
5105  _bnd_mat_side_cache[tid][bnd_id] = false;
5106 
5107  if (_nl->needMaterialOnSide(bnd_id, tid) || _aux->needMaterialOnSide(bnd_id))
5108  _bnd_mat_side_cache[tid][bnd_id] = true;
5109 
5110  else if (_side_user_objects.hasActiveBoundaryObjects(bnd_id, tid))
5111  _bnd_mat_side_cache[tid][bnd_id] = true;
5112  }
5113 
5114  return _bnd_mat_side_cache[tid][bnd_id];
5115 }
5116 
5117 bool
5119 {
5120  if (_block_mat_side_cache[tid].find(subdomain_id) == _block_mat_side_cache[tid].end())
5121  {
5122  _block_mat_side_cache[tid][subdomain_id] = false;
5123 
5124  if (_nl->needMaterialOnSide(subdomain_id, tid))
5125  _block_mat_side_cache[tid][subdomain_id] = true;
5126  else if (_internal_side_user_objects.hasActiveBlockObjects(subdomain_id, tid))
5127  _block_mat_side_cache[tid][subdomain_id] = true;
5128  }
5129 
5130  return _block_mat_side_cache[tid][subdomain_id];
5131 }
5132 
5133 bool
5135 {
5136  return _needs_old_newton_iter;
5137 }
5138 
5139 void
5141 {
5142  _needs_old_newton_iter = state;
5143 }
const std::string & name() const
Get the name of the object.
Definition: MooseObject.h:47
Object is evaluated in every residual computation.
Definition: MooseTypes.h:96
bool _reinit_displaced_elem
This class determines the maximum number of Quadrature Points and Shape Functions used for a given si...
Definition: MaxQpsThread.h:33
virtual void addResidual(THREAD_ID tid) override
std::string indent(unsigned int spaces)
Create empty string for indenting.
Definition: ConsoleUtils.C:34
Interface for objects that need parallel consistent random numbers without patterns over the course o...
VectorPostprocessorValue & declareVector(const std::string &vpp_name, const std::string &vector_name)
Initialization method, sets the current and old value to 0.0 for this VectorPostprocessor.
unsigned int max() const
Definition: MaxQpsThread.h:45
bool adaptMesh(std::string marker_name=std::string())
Adapts the mesh based on the error estimator used.
Definition: Adaptivity.C:125
void sort(THREAD_ID tid=0)
Sort the objects using the DependencyResolver.
bool initialAdaptMesh()
Used during initial adaptivity.
Definition: Adaptivity.C:199
void addMultiApp(const std::string &multi_app_name, const std::string &name, InputParameters parameters)
Add a MultiApp to the problem.
const std::vector< dof_id_type > & allDofIndices() const
Get all global dofindices for the variable.
virtual void setActiveElementalMooseVariables(const std::set< MooseVariable * > &moose_vars, THREAD_ID tid) override
Set the MOOSE variables to be reinited on each element.
MeshBase::const_element_iterator activeLocalElementsBegin()
Calls active_local_nodes_begin/end() on the underlying libMesh mesh object.
Definition: MooseMesh.C:2065
RealVectorValue RealGradient
Definition: Assembly.h:43
Helper class for holding the preconditioning blocks to fill.
virtual const std::set< unsigned int > & getActiveMaterialProperties(THREAD_ID tid) override
Get the material properties required by the current computing thread.
virtual void setResidualNeighbor(NumericVector< Number > &residual, THREAD_ID tid) override
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...
void outputStep(ExecFlagType type)
Calls the outputStep method for each output object.
std::shared_ptr< T > getActiveObject(const std::string &name, THREAD_ID tid=0) const
ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:724
bool _requires_nonlocal_coupling
nonlocal coupling requirement flag
Definition: SubProblem.h:461
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:46
virtual void prepare(const Elem *elem, THREAD_ID tid) override
void shift()
Shift the material properties in time.
Base class for function objects.
Definition: Function.h:46
const std::vector< std::pair< std::string, VectorPostprocessorState > > & vectors(const std::string &vpp_name) const
Get the map of vectors for a particular VectorPostprocessor.
void initPetscOutput()
Reinitialize petsc output for proper linear/nonlinear iteration display.
virtual System & getSystem(const std::string &var_name) override
Returns the equation system containing the variable provided.
void mooseWarning(Args &&...args) const
Definition: MooseObject.h:89
std::shared_ptr< Postprocessor > getPostprocessorPointer(std::shared_ptr< MooseObject > mo)
Small helper function used by addPostprocessor to try to get a Postprocessor pointer from a MooseObje...
virtual void addResidualNeighbor(THREAD_ID tid) override
A class for creating restricted objects.
Definition: Restartable.h:31
virtual void computeNullSpace(NonlinearImplicitSystem &sys, std::vector< NumericVector< Number > * > &sp)
std::vector< std::shared_ptr< MaterialData > > _neighbor_material_data
Factory & _factory
The Factory for building objects.
Definition: SubProblem.h:418
unsigned int _rz_coord_axis
Storage for RZ axis selection.
Definition: SubProblem.h:467
MaterialPropertyStorage & _bnd_material_props
bool hasUserObject(const std::string &name)
Check if there if a user object of given name.
virtual void addGhostedElem(dof_id_type elem_id) override
Will make sure that all dofs connected to elem_id are ghosted to this processor.
void setNonlocalCouplingMatrix()
Set custom coupling matrix for variables requiring nonlocal contribution.
Threads::spin_mutex get_function_mutex
unsigned int getInitialSteps() const
Pull out the number of initial steps previously set by calling init()
Definition: Adaptivity.h:94
void initialSetup()
Calls the initialSetup function for each of the output objects.
SolverParams _solver_params
void execMultiAppTransfers(ExecFlagType type, MultiAppTransfer::DIRECTION direction)
Execute MultiAppTransfers associate with execution flag and direction.
SolverParams & solverParams()
Get the solver parameters.
virtual void prepareFace(const Elem *elem, THREAD_ID tid) override
virtual void setInputParametersFEProblem(InputParameters &parameters)
ExecFlagType _current_execute_on_flag
Current execute_on flag.
virtual void addPredictor(const std::string &type, const std::string &name, InputParameters parameters)
virtual void prepareMaterials(SubdomainID blk_id, THREAD_ID tid)
Add the MooseVariables that the current materials depend on to the dependency list.
bool duplicateVariableCheck(const std::string &var_name, const FEType &type, bool is_aux)
Helper to check for duplicate variable names across systems or within a single system.
Keeps track of stuff related to assembling.
Definition: Assembly.h:63
void addMaterial(const std::string &kernel_name, const std::string &name, InputParameters parameters)
void storeValue(const std::string &name, PostprocessorValue value)
virtual void reinitElemPhys(const Elem *elem, std::vector< Point > phys_points_in_elem, THREAD_ID tid) override
void updateBlockMatPropDependency(SubdomainID id, std::set< unsigned int > &needed_mat_props, THREAD_ID tid=0) const
Class for stuff related to variables.
Definition: MooseVariable.h:43
bool & setFileRestart()
This method is here so we can determine whether or not we need to use a separate reader to read the m...
Definition: MooseApp.h:271
NonlinearSystemBase & getNonlinearSystemBase()
virtual void addVectorPostprocessor(std::string pp_name, const std::string &name, InputParameters parameters)
std::map< SubdomainID, Moose::CoordinateSystemType > _coord_sys
nonlocal coupling matrix;
Definition: SubProblem.h:423
void addScalarVariable(const std::string &var_name, Order order, Real scale_factor=1., const std::set< SubdomainID > *const active_subdomains=NULL)
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:2115
bool isRecovering() const
Whether or not this is a "recover" calculation.
Definition: MooseApp.C:607
std::shared_ptr< NonlinearSystemBase > _nl
bool _has_jacobian
Indicates if the Jacobian was computed.
void setSolverDefaults(FEProblemBase &problem)
Definition: Moose.C:1233
void addPrivateParam(const std::string &name, const T &value)
These method add a parameter to the InputParameters object which can be retrieved like any other para...
std::pair< typename M::iterator, bool > moose_try_emplace(M &m, const typename M::key_type &k, Args &&...args)
Function to mirror the behavior of the C++17 std::map::try_emplace() method (no hint).
Definition: Moose.h:66
bool _has_dampers
Whether or not this system has any Dampers associated with it.
virtual Distribution & getDistribution(const std::string &name)
unsigned int multiAppLevel() const
The MultiApp Level.
Definition: MooseApp.h:465
void computeTransientImplicitResidual(Real time, const NumericVector< Number > &u, const NumericVector< Number > &udot, NumericVector< Number > &residual)
Evaluates transient residual G in canonical semidiscrete form G(t,U,Udot) = F(t,U) ...
bool _has_nonlocal_coupling
Indicates if nonlocal coupling is required/exists.
virtual void reinitNodeNeighbor(const Node *node, THREAD_ID tid) override
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:490
This is the base class for Samplers.
Definition: Sampler.h:46
virtual bool reinitDirac(const Elem *elem, THREAD_ID tid) override
Returns true if the Problem has Dirac kernels it needs to compute on elem.
Object is evaluated only once at the beginning of the simulation.
Definition: MooseTypes.h:94
virtual void computeResidualType(const NumericVector< Number > &soln, NumericVector< Number > &residual, Moose::KernelType type=Moose::KT_ALL)
Base class for predictors.
Definition: Predictor.h:39
void addDiracKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
subdomain_id_type SubdomainID
Definition: MooseTypes.h:77
virtual void postExecute()
Method called at the end of the simulation.
SubdomainID getSubdomainID(const SubdomainName &subdomain_name) const
Get the associated subdomain ID for the subdomain name.
Definition: MooseMesh.C:1036
void petscSetDefaults(FEProblemBase &problem)
Sets the default options for PETSc.
Definition: PetscSupport.C:560
InputParameters getValidParams(const std::string &name)
Get valid parameters for the object.
Definition: Factory.C:26
InputParameterWarehouse & getInputParameterWarehouse()
Get the InputParameterWarehouse for MooseObjects.
Definition: MooseApp.C:1121
void addAuxKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
ExecuteMooseObjectWarehouse< Control > _control_warehouse
The control logic warehouse.
AuxGroupExecuteMooseObjectWarehouse< NodalUserObject > _nodal_user_objects
void setCoupling(Moose::CouplingType type)
Set the coupling between variables TODO: allow user-defined coupling.
unsigned int _cycles_completed
virtual void init() override
const ExecFlagType & getCurrentExecuteOnFlag() const
Return the current execution flag.
virtual void prepareNeighborShapes(unsigned int var, THREAD_ID tid) override
const std::map< BoundaryID, std::vector< std::shared_ptr< T > > > & getActiveBoundaryObjects(THREAD_ID tid=0) const
void parentOutputPositionChanged()
Calls parentOutputPositionChanged() on all sub apps.
const std::set< SubdomainID > & meshSubdomains() const
Returns a read-only reference to the set of subdomains currently present in the Mesh.
Definition: MooseMesh.C:2157
PostprocessorValue & getPostprocessorValueOld(const std::string &name)
Get the reference to the old value of a post-processor.
VectorPostprocessorValue & getVectorPostprocessorValueOld(const VectorPostprocessorName &vpp_name, const std::string &vector_name)
The the old value of an post-processor.
std::string getRecoverFileBase()
The file_base for the recovery file.
Definition: MooseApp.h:318
NonlocalIntegratedBC is used for solving integral terms in integro-differential equations.
bool hasActiveObject(const std::string &name, THREAD_ID tid=0) const
Convenience functions for checking/getting specific objects.
AuxGroupExecuteMooseObjectWarehouse< ElementUserObject > _elemental_user_objects
void serializeSolution()
virtual void reinitMaterialsNeighbor(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
void updateErrorVectors()
Update the ErrorVectors that have been requested through calls to getErrorVector().
Definition: Adaptivity.C:285
virtual void initNullSpaceVectors(const InputParameters &parameters, NonlinearSystemBase &nl)
virtual void computeMarkers()
void reportMooseObjectDependency(MooseObject *a, MooseObject *b)
Register a MOOSE object dependency so we can either order operations properly or report when we canno...
void indentMessage(const std::string &prefix, std::string &message, const char *color=COLOR_CYAN)
Indents the supplied message given the prefix and color.
Definition: MooseUtils.C:438
MaterialDataType
MaterialData types.
Definition: MooseTypes.h:129
void cacheChangedLists()
Cache information about what elements were refined and coarsened in the previous step.
Definition: MooseMesh.C:520
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
void updateActiveObjects()
Update the active objects in the warehouses.
Stores the stateful material properties computed by materials.
virtual void computeIndicatorsAndMarkers()
Definition: Marker.h:42
bool getRecomputeMarkersFlag() const
Pull out the _recompute_markers_during_cycles flag previously set through the AdaptivityAction.
Definition: Adaptivity.h:121
virtual Real & time() const
virtual void setException(const std::string &message)
Set an exception.
void sort(THREAD_ID tid=0)
Performs a sort using the DependencyResolver.
virtual void jacobianSetup(THREAD_ID tid=0) const
virtual bool hasActiveMaterialProperties(THREAD_ID tid)
Method to check whether or not a list of active material roperties has been set.
Definition: SubProblem.C:93
std::shared_ptr< Material > getMaterial(std::string name, Moose::MaterialDataType type, THREAD_ID tid=0, bool no_warn=false)
Return a pointer to a Material object.
bool hasVectorPostprocessor(const std::string &name)
Returns a true value if the VectorPostprocessor exists.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
bool hasPostprocessor(const std::string &name)
Check existence of the postprocessor.
void init(const std::string &name)
Initialization method, sets the current and old value to 0.0 for this postprocessor.
virtual void computeBounds(NonlinearImplicitSystem &sys, NumericVector< Number > &lower, NumericVector< Number > &upper)
RealTensorValue RealTensor
Definition: Assembly.h:46
MultiApp Implementation for Transient Apps.
void addNodalKernel(const std::string &kernel_name, const std::string &name, InputParameters parameters)
virtual void getDiracElements(std::set< const Elem * > &elems) override
Fills "elems" with the elements that should be looped over for Dirac Kernels.
bool needMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid)
These methods are used to determine whether stateful material properties need to be stored on interna...
std::vector< VariableSecond > _second_zero
virtual const std::set< unsigned int > & getActiveMaterialProperties(THREAD_ID tid)
Get the material properties required by the current computing thread.
Definition: SubProblem.C:87
virtual void onTimestepEnd() override
All Distributions should inherit from this class.
Definition: Distribution.h:28
void setCouplingMatrix(std::unique_ptr< CouplingMatrix > cm)
Set custom coupling matrix.
virtual bool updateSolution(NumericVector< Number > &vec_solution, NumericVector< Number > &ghosted_solution)
Update the solution.
bool _has_exception
Whether or not an exception has occurred.
bool _needs_old_newton_iter
Indicates that we need to compute variable values for previous Newton iteration.
bool haveXFEM()
Find out whether the current analysis is using XFEM.
void registerRandomInterface(RandomInterface &random_interface, const std::string &name)
bool _has_time_integrator
Indicates whether or not this executioner has a time integrator (during setup)
virtual bool hasActiveElementalMooseVariables(THREAD_ID tid)
Whether or not a list of active elemental moose variables has been set.
Definition: SubProblem.C:67
bool hasVectorPostprocessor(const std::string &name)
Check existence of the VectorPostprocessor.
ExecuteMooseObjectWarehouse< TransientMultiApp > _transient_multi_apps
Storage for TransientMultiApps (only needed for calling &#39;computeDT&#39;)
virtual void addObject(std::shared_ptr< T > object, THREAD_ID tid=0)
Adds an object to the storage structure.
void errorIfDistributedMesh(std::string name) const
Generate a unified error message if the underlying libMesh mesh is a DistributedMesh.
Definition: MooseMesh.C:2526
VectorPostprocessorData & getVectorPostprocessorData()
This class is here to combine the Postprocessor interface and the base class Postprocessor object alo...
ExecuteMooseObjectWarehouse< Transfer > _from_multi_app_transfers
Transfers executed just after MultiApps to transfer data from them.
void setVariableAllDoFMap(const std::vector< MooseVariable * > moose_vars)
Base class for a system (of equations)
Definition: SystemBase.h:91
const std::map< unsigned int, std::string > statefulPropNames() const
PostprocessorValue & getPostprocessorValueOlder(const std::string &name)
Get the reference to the older value of a post-processor.
virtual void reinitElem(const Elem *elem, THREAD_ID tid) override
virtual void reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
Definition: