www.mooseframework.org
FEProblemBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "FEProblemBase.h"
11 #include "AuxiliarySystem.h"
13 #include "MooseEnum.h"
14 #include "Resurrector.h"
15 #include "Factory.h"
16 #include "MooseUtils.h"
17 #include "DisplacedProblem.h"
18 #include "SystemBase.h"
19 #include "MaterialData.h"
25 #include "ComputeIndicatorThread.h"
26 #include "ComputeMarkerThread.h"
29 #include "MaxQpsThread.h"
30 #include "ActionWarehouse.h"
31 #include "Conversion.h"
32 #include "Material.h"
33 #include "ConstantIC.h"
34 #include "Parser.h"
35 #include "ElementH1Error.h"
36 #include "Function.h"
37 #include "NonlinearSystem.h"
38 #include "Distribution.h"
39 #include "Sampler.h"
40 #include "PetscSupport.h"
41 #include "RandomInterface.h"
42 #include "RandomData.h"
43 #include "MooseEigenSystem.h"
44 #include "MooseParsedFunction.h"
45 #include "MeshChangedInterface.h"
47 #include "ScalarInitialCondition.h"
48 #include "ElementPostprocessor.h"
49 #include "NodalPostprocessor.h"
50 #include "SidePostprocessor.h"
52 #include "InterfacePostprocessor.h"
53 #include "GeneralPostprocessor.h"
59 #include "Indicator.h"
60 #include "Marker.h"
61 #include "MultiApp.h"
62 #include "MultiAppTransfer.h"
63 #include "TransientMultiApp.h"
64 #include "ElementUserObject.h"
65 #include "NodalUserObject.h"
66 #include "SideUserObject.h"
67 #include "InternalSideUserObject.h"
68 #include "InterfaceUserObject.h"
69 #include "GeneralUserObject.h"
71 #include "InternalSideIndicator.h"
72 #include "Transfer.h"
73 #include "MultiAppTransfer.h"
74 #include "MultiMooseEnum.h"
75 #include "Predictor.h"
76 #include "Assembly.h"
77 #include "Control.h"
78 #include "XFEMInterface.h"
79 #include "ConsoleUtils.h"
80 #include "NonlocalKernel.h"
81 #include "NonlocalIntegratedBC.h"
82 #include "ShapeElementUserObject.h"
83 #include "ShapeSideUserObject.h"
84 #include "MooseVariableFE.h"
85 #include "MooseVariableScalar.h"
87 #include "TimeIntegrator.h"
88 #include "LineSearch.h"
90 
91 #include "libmesh/exodusII_io.h"
92 #include "libmesh/quadrature.h"
93 #include "libmesh/coupling_matrix.h"
94 #include "libmesh/nonlinear_solver.h"
95 #include "libmesh/sparse_matrix.h"
96 
97 // Anonymous namespace for helper function
98 namespace
99 {
103 bool
104 sortMooseVariables(MooseVariableFEBase * a, MooseVariableFEBase * b)
105 {
106  return a->number() < b->number();
107 }
108 } // namespace
109 
110 Threads::spin_mutex get_function_mutex;
111 
112 template <>
115 {
117  params.addParam<unsigned int>("null_space_dimension", 0, "The dimension of the nullspace");
118  params.addParam<unsigned int>(
119  "transpose_null_space_dimension", 0, "The dimension of the transpose nullspace");
120  params.addParam<unsigned int>(
121  "near_null_space_dimension", 0, "The dimension of the near nullspace");
122  params.addParam<bool>("solve",
123  true,
124  "Whether or not to actually solve the Nonlinear system. "
125  "This is handy in the case that all you want to do is "
126  "execute AuxKernels, Transfers, etc. without actually "
127  "solving anything");
128  params.addParam<bool>("use_nonlinear",
129  true,
130  "Determines whether to use a Nonlinear vs a "
131  "Eigenvalue system (Automatically determined based "
132  "on executioner)");
133  params.addParam<bool>("error_on_jacobian_nonzero_reallocation",
134  false,
135  "This causes PETSc to error if it had to reallocate memory in the Jacobian "
136  "matrix due to not having enough nonzeros");
137  params.addParam<bool>("ignore_zeros_in_jacobian",
138  false,
139  "Do not explicitly store zero values in "
140  "the Jacobian matrix if true");
141  params.addParam<bool>("force_restart",
142  false,
143  "EXPERIMENTAL: If true, a sub_app may use a "
144  "restart file instead of using of using the master "
145  "backup file");
146  params.addParam<bool>("skip_additional_restart_data",
147  false,
148  "True to skip additional data in equation system for restart. It is useful "
149  "for starting a transient calculation with a steady-state solution");
150  params.addParam<bool>("skip_nl_system_check",
151  false,
152  "True to skip the NonlinearSystem check for work to do (e.g. Make sure "
153  "that there are variables to solve for).");
154 
160  params.addParam<std::vector<SubdomainName>>("block", "Block IDs for the coordinate systems");
161  MultiMooseEnum coord_types("XYZ RZ RSPHERICAL", "XYZ");
162  MooseEnum rz_coord_axis("X=0 Y=1", "Y");
163  params.addParam<MultiMooseEnum>(
164  "coord_type", coord_types, "Type of the coordinate system per block param");
165  params.addParam<MooseEnum>(
166  "rz_coord_axis", rz_coord_axis, "The rotation axis (X | Y) for axisymetric coordinates");
167  params.addParam<bool>(
168  "kernel_coverage_check", true, "Set to false to disable kernel->subdomain coverage check");
169  params.addParam<bool>("material_coverage_check",
170  true,
171  "Set to false to disable material->subdomain coverage check");
172  params.addParam<bool>("parallel_barrier_messaging",
173  false,
174  "Displays messaging from parallel "
175  "barrier notifications when executing "
176  "or transferring to/from Multiapps "
177  "(default: false)");
178 
179  params.addParam<FileNameNoExtension>("restart_file_base",
180  "File base name used for restart (e.g. "
181  "<path>/<filebase> or <path>/LATEST to "
182  "grab the latest file available)");
183 
184  params.addParam<std::vector<TagName>>("extra_tag_vectors",
185  "Extra vectors to add to the system that can be filled by "
186  "objects which compute residuals and Jacobians (Kernels, "
187  "BCs, etc.) by setting tags on them.");
188 
189  params.addParam<std::vector<TagName>>("extra_tag_matrices",
190  "Extra matrices to add to the system that can be filled "
191  "by objects which compute residuals and Jacobians "
192  "(Kernels, BCs, etc.) by setting tags on them.");
193 
194  return params;
195 }
196 
198  : SubProblem(parameters),
199  Restartable(this, "FEProblemBase"),
200  _mesh(*getCheckedPointerParam<MooseMesh *>("mesh")),
201  _eq(_mesh),
202  _initialized(false),
203  _solve(getParam<bool>("solve")),
204  _transient(false),
205  _time(declareRestartableData<Real>("time")),
206  _time_old(declareRestartableData<Real>("time_old")),
207  _t_step(declareRecoverableData<int>("t_step")),
208  _dt(declareRestartableData<Real>("dt")),
209  _dt_old(declareRestartableData<Real>("dt_old")),
210  _nl(nullptr),
211  _aux(nullptr),
212  _coupling(Moose::COUPLING_DIAG),
213  _distributions(/*threaded=*/false),
214  _samplers(_app.getExecuteOnEnum()),
215  _material_props(
216  declareRestartableDataWithContext<MaterialPropertyStorage>("material_props", &_mesh)),
217  _bnd_material_props(
218  declareRestartableDataWithContext<MaterialPropertyStorage>("bnd_material_props", &_mesh)),
219  _neighbor_material_props(declareRestartableDataWithContext<MaterialPropertyStorage>(
220  "neighbor_material_props", &_mesh)),
221  _pps_data(*this),
222  _vpps_data(*this),
223  // TODO: delete the following line after apps have been updated to not call getUserObjects
224  _all_user_objects(_app.getExecuteOnEnum()),
225  _multi_apps(_app.getExecuteOnEnum()),
226  _transient_multi_apps(_app.getExecuteOnEnum()),
227  _transfers(_app.getExecuteOnEnum(), /*threaded=*/false),
228  _to_multi_app_transfers(_app.getExecuteOnEnum(), /*threaded=*/false),
229  _from_multi_app_transfers(_app.getExecuteOnEnum(), /*threaded=*/false),
230 #ifdef LIBMESH_ENABLE_AMR
231  _adaptivity(*this),
232  _cycles_completed(0),
233 #endif
234  _displaced_mesh(nullptr),
235  _geometric_search_data(*this, _mesh),
236  _mortar_data(),
237  _reinit_displaced_elem(false),
238  _reinit_displaced_face(false),
239  _input_file_saved(false),
240  _has_dampers(false),
241  _has_constraints(false),
242  _snesmf_reuse_base(true),
243  _snesmf_reuse_base_set_by_user(false),
244  _has_initialized_stateful(false),
245  _const_jacobian(false),
246  _has_jacobian(false),
247  _needs_old_newton_iter(false),
248  _has_nonlocal_coupling(false),
249  _calculate_jacobian_in_uo(false),
250  _kernel_coverage_check(getParam<bool>("kernel_coverage_check")),
251  _material_coverage_check(getParam<bool>("material_coverage_check")),
252  _max_qps(std::numeric_limits<unsigned int>::max()),
253  _max_shape_funcs(std::numeric_limits<unsigned int>::max()),
254  _max_scalar_order(INVALID_ORDER),
255  _has_time_integrator(false),
256  _has_exception(false),
257  _parallel_barrier_messaging(getParam<bool>("parallel_barrier_messaging")),
258  _current_execute_on_flag(EXEC_NONE),
259  _control_warehouse(_app.getExecuteOnEnum(), /*threaded=*/false),
260  _line_search(nullptr),
261  _using_ad_mat_props(false),
262  _error_on_jacobian_nonzero_reallocation(
263  getParam<bool>("error_on_jacobian_nonzero_reallocation")),
264  _ignore_zeros_in_jacobian(getParam<bool>("ignore_zeros_in_jacobian")),
265  _force_restart(getParam<bool>("force_restart")),
266  _skip_additional_restart_data(getParam<bool>("skip_additional_restart_data")),
267  _skip_nl_system_check(getParam<bool>("skip_nl_system_check")),
268  _fail_next_linear_convergence_check(false),
269  _started_initial_setup(false),
270  _has_internal_edge_residual_objects(false),
271  _initial_setup_timer(registerTimedSection("initialSetup", 2)),
272  _project_solution_timer(registerTimedSection("projectSolution", 2)),
273  _compute_indicators_timer(registerTimedSection("computeIndicators", 1)),
274  _compute_markers_timer(registerTimedSection("computeMarkers", 1)),
275  _compute_user_objects_timer(registerTimedSection("computeUserObjects", 1)),
276  _execute_controls_timer(registerTimedSection("executeControls", 1)),
277  _execute_samplers_timer(registerTimedSection("executeSamplers", 1)),
278  _update_active_objects_timer(registerTimedSection("updateActiveObjects", 5)),
279  _reinit_because_of_ghosting_or_new_geom_objects_timer(
280  registerTimedSection("reinitBecauseOfGhostingOrNewGeomObjects", 3)),
281  _exec_multi_app_transfers_timer(registerTimedSection("execMultiAppTransfers", 1)),
282  _init_timer(registerTimedSection("init", 2)),
283  _eq_init_timer(registerTimedSection("EquationSystems::Init", 2)),
284  _solve_timer(registerTimedSection("solve", 1)),
285  _check_exception_and_stop_solve_timer(registerTimedSection("checkExceptionAndStopSolve", 5)),
286  _advance_state_timer(registerTimedSection("advanceState", 5)),
287  _restore_solutions_timer(registerTimedSection("restoreSolutions", 5)),
288  _save_old_solutions_timer(registerTimedSection("saveOldSolutions", 5)),
289  _restore_old_solutions_timer(registerTimedSection("restoreOldSolutions", 5)),
290  _output_step_timer(registerTimedSection("outputStep", 1)),
291  _on_timestep_begin_timer(registerTimedSection("onTimestepBegin", 2)),
292  _compute_residual_l2_norm_timer(registerTimedSection("computeResidualL2Norm", 2)),
293  _compute_residual_sys_timer(registerTimedSection("computeResidualSys", 5)),
294  _compute_residual_internal_timer(registerTimedSection("computeResidualInternal", 1)),
295  _compute_residual_type_timer(registerTimedSection("computeResidualType", 5)),
296  _compute_transient_implicit_residual_timer(
297  registerTimedSection("computeTransientImplicitResidual", 2)),
298  _compute_residual_tags_timer(registerTimedSection("computeResidualTags", 5)),
299  _compute_jacobian_internal_timer(registerTimedSection("computeJacobianInternal", 1)),
300  _compute_jacobian_tags_timer(registerTimedSection("computeJacobianTags", 5)),
301  _compute_jacobian_blocks_timer(registerTimedSection("computeTransientImplicitJacobian", 2)),
302  _compute_bounds_timer(registerTimedSection("computeBounds", 1)),
303  _compute_post_check_timer(registerTimedSection("computePostCheck", 2)),
304  _compute_damping_timer(registerTimedSection("computeDamping", 1)),
305  _possibly_rebuild_geom_search_patches_timer(
306  registerTimedSection("possiblyRebuildGeomSearchPatches", 5)),
307  _initial_adapt_mesh_timer(registerTimedSection("initialAdaptMesh", 2)),
308  _adapt_mesh_timer(registerTimedSection("adaptMesh", 3)),
309  _update_mesh_xfem_timer(registerTimedSection("updateMeshXFEM", 5)),
310  _mesh_changed_timer(registerTimedSection("meshChanged", 3)),
311  _mesh_changed_helper_timer(registerTimedSection("meshChangedHelper", 5)),
312  _check_problem_integrity_timer(registerTimedSection("notifyWhenMeshChanges", 5)),
313  _serialize_solution_timer(registerTimedSection("serializeSolution", 3)),
314  _check_nonlinear_convergence_timer(registerTimedSection("checkNonlinearConvergence", 5)),
315  _check_linear_convergence_timer(registerTimedSection("checkLinearConvergence", 5)),
316  _update_geometric_search_timer(registerTimedSection("updateGeometricSearch", 3)),
317  _exec_multi_apps_timer(registerTimedSection("execMultiApps", 1)),
318  _backup_multi_apps_timer(registerTimedSection("backupMultiApps", 5)),
319  _u_dot_requested(false),
320  _u_dotdot_requested(false),
321  _u_dot_old_requested(false),
322  _u_dotdot_old_requested(false)
323 {
324  // Possibly turn off default ghosting in libMesh
325  _eq.enable_default_ghosting(_default_ghosting);
326 
327  _time = 0.0;
328  _time_old = 0.0;
329  _t_step = 0;
330  _dt = 0;
331  _dt_old = _dt;
332 
333  unsigned int n_threads = libMesh::n_threads();
334 
335  _real_zero.resize(n_threads, 0.);
336  _scalar_zero.resize(n_threads);
337  _zero.resize(n_threads);
338  _ad_zero.resize(n_threads);
339  _grad_zero.resize(n_threads);
340  _ad_grad_zero.resize(n_threads);
341  _second_zero.resize(n_threads);
342  _ad_second_zero.resize(n_threads);
343  _second_phi_zero.resize(n_threads);
344  _point_zero.resize(n_threads);
345  _vector_zero.resize(n_threads);
346  _vector_curl_zero.resize(n_threads);
347  _uo_jacobian_moose_vars.resize(n_threads);
348 
349  _material_data.resize(n_threads);
350  _bnd_material_data.resize(n_threads);
351  _neighbor_material_data.resize(n_threads);
352  for (unsigned int i = 0; i < n_threads; i++)
353  {
354  _material_data[i] = std::make_shared<MaterialData>(_material_props);
355  _bnd_material_data[i] = std::make_shared<MaterialData>(_bnd_material_props);
356  _neighbor_material_data[i] = std::make_shared<MaterialData>(_neighbor_material_props);
357  }
358 
359  _active_elemental_moose_variables.resize(n_threads);
360 
361  _block_mat_side_cache.resize(n_threads);
362  _bnd_mat_side_cache.resize(n_threads);
363 
364  _resurrector = libmesh_make_unique<Resurrector>(*this);
365 
366  _eq.parameters.set<FEProblemBase *>("_fe_problem_base") = this;
367 
368  setCoordSystem(getParam<std::vector<SubdomainName>>("block"),
369  getParam<MultiMooseEnum>("coord_type"));
370  setAxisymmetricCoordAxis(getParam<MooseEnum>("rz_coord_axis"));
371 
372  if (isParamValid("restart_file_base"))
373  {
374  std::string restart_file_base = getParam<FileNameNoExtension>("restart_file_base");
375  restart_file_base = MooseUtils::convertLatestCheckpoint(restart_file_base);
376  _console << "\nUsing " << restart_file_base << " for restart.\n\n";
377  setRestartFile(restart_file_base);
378  }
379 }
380 
381 void
383 {
384  // add vectors and their tags to system
385  auto & vectors = getParam<std::vector<TagName>>("extra_tag_vectors");
386  for (auto & vector : vectors)
387  {
388  auto tag = addVectorTag(vector);
389  _nl->addVector(tag, false, GHOSTED);
390  }
391 
392  // add matrices and their tags
393  auto & matrices = getParam<std::vector<TagName>>("extra_tag_matrices");
394  for (auto & matrix : matrices)
395  {
396  auto tag = addMatrixTag(matrix);
397  _nl->addMatrix(tag);
398  }
399 }
400 
401 void
403 {
404  unsigned int n_threads = libMesh::n_threads();
405 
406  _assembly.resize(n_threads);
407  for (unsigned int i = 0; i < n_threads; ++i)
408  _assembly[i] = libmesh_make_unique<Assembly>(nl, i);
409 }
410 
411 void
413 {
414  unsigned int dimNullSpace = parameters.get<unsigned int>("null_space_dimension");
415  unsigned int dimTransposeNullSpace =
416  parameters.get<unsigned int>("transpose_null_space_dimension");
417  unsigned int dimNearNullSpace = parameters.get<unsigned int>("near_null_space_dimension");
418  for (unsigned int i = 0; i < dimNullSpace; ++i)
419  {
420  std::ostringstream oss;
421  oss << "_" << i;
422  // do not project, since this will be recomputed, but make it ghosted, since the near nullspace
423  // builder might march over all nodes
424  nl.addVector("NullSpace" + oss.str(), false, GHOSTED);
425  }
426  _subspace_dim["NullSpace"] = dimNullSpace;
427  for (unsigned int i = 0; i < dimTransposeNullSpace; ++i)
428  {
429  std::ostringstream oss;
430  oss << "_" << i;
431  // do not project, since this will be recomputed, but make it ghosted, since the near nullspace
432  // builder might march over all nodes
433  nl.addVector("TransposeNullSpace" + oss.str(), false, GHOSTED);
434  }
435  _subspace_dim["TransposeNullSpace"] = dimTransposeNullSpace;
436  for (unsigned int i = 0; i < dimNearNullSpace; ++i)
437  {
438  std::ostringstream oss;
439  oss << "_" << i;
440  // do not project, since this will be recomputed, but make it ghosted, since the near-nullspace
441  // builder might march over all semilocal nodes
442  nl.addVector("NearNullSpace" + oss.str(), false, GHOSTED);
443  }
444  _subspace_dim["NearNullSpace"] = dimNearNullSpace;
445 }
446 
448 {
449  // Flush the Console stream, the underlying call to Console::mooseConsole
450  // relies on a call to Output::checkInterval that has references to
451  // _time, etc. If it is not flushed here memory problems arise if you have
452  // an unflushed stream and start destructing things.
453  _console << std::flush;
454 
455  unsigned int n_threads = libMesh::n_threads();
456  for (unsigned int i = 0; i < n_threads; i++)
457  {
458  _zero[i].release();
459  _scalar_zero[i].release();
460  _grad_zero[i].release();
461  _second_zero[i].release();
462  _second_phi_zero[i].release();
463  _vector_zero[i].release();
464  _vector_curl_zero[i].release();
465  _ad_zero[i].release();
466  _ad_grad_zero[i].release();
467  _ad_second_zero[i].release();
468  }
469 }
470 
473 {
474  std::map<SubdomainID, Moose::CoordinateSystemType>::iterator it = _coord_sys.find(sid);
475  if (it != _coord_sys.end())
476  return (*it).second;
477  else
478  mooseError("Requested subdomain ", sid, " does not exist.");
479 }
480 
481 void
482 FEProblemBase::setCoordSystem(const std::vector<SubdomainName> & blocks,
483  const MultiMooseEnum & coord_sys)
484 {
485  const std::set<SubdomainID> & subdomains = _mesh.meshSubdomains();
486  if (blocks.size() == 0)
487  {
488  // no blocks specified -> assume the whole domain
489  Moose::CoordinateSystemType coord_type = Moose::COORD_XYZ; // all is going to be XYZ by default
490  if (coord_sys.size() == 0)
491  ; // relax, do nothing
492  else if (coord_sys.size() == 1)
493  coord_type = Moose::stringToEnum<Moose::CoordinateSystemType>(
494  coord_sys[0]); // one system specified, the whole domain is going to have that system
495  else
496  mooseError("Multiple coordinate systems specified, but no blocks given.");
497 
498  for (const auto & sbd : subdomains)
499  _coord_sys[sbd] = coord_type;
500  }
501  else
502  {
503  // user specified 'blocks' but not coordinate systems
504  if (coord_sys.size() == 0)
505  {
506  // set all blocks to cartesian coordinate system
507  for (const auto & block : blocks)
508  {
509  SubdomainID sid = _mesh.getSubdomainID(block);
511  }
512  }
513  else if (coord_sys.size() == 1)
514  {
515  // set all blocks to the coordinate system specified by `coord_sys[0]`
516  Moose::CoordinateSystemType coord_type =
517  Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[0]);
518  for (const auto & block : blocks)
519  {
520  SubdomainID sid = _mesh.getSubdomainID(block);
521  _coord_sys[sid] = coord_type;
522  }
523  }
524  else
525  {
526  if (blocks.size() != coord_sys.size())
527  mooseError("Number of blocks and coordinate systems does not match.");
528 
529  for (unsigned int i = 0; i < blocks.size(); i++)
530  {
531  SubdomainID sid = _mesh.getSubdomainID(blocks[i]);
532  Moose::CoordinateSystemType coord_type =
533  Moose::stringToEnum<Moose::CoordinateSystemType>(coord_sys[i]);
534  _coord_sys[sid] = coord_type;
535  }
536 
537  for (const auto & sid : subdomains)
538  if (_coord_sys.find(sid) == _coord_sys.end())
539  mooseError("Subdomain '" + Moose::stringify(sid) +
540  "' does not have a coordinate system specified.");
541  }
542  }
543 }
544 
545 void
547 {
548  _rz_coord_axis = rz_coord_axis;
549 }
550 
551 void
553 {
554  _nl->addExtraVectors();
555  _aux->addExtraVectors();
556 }
557 
558 const ConstElemRange &
560 {
562  {
564  libmesh_make_unique<ConstElemRange>(_mesh.getMesh().evaluable_elements_begin(_nl->dofMap()),
565  _mesh.getMesh().evaluable_elements_end(_nl->dofMap()));
566  }
568 }
569 
570 void
572 {
573  TIME_SECTION(_initial_setup_timer);
574 
575  // set state flag indicating that we are in or beyond initialSetup.
576  // This can be used to throw errors in methods that _must_ be called at construction time.
577  _started_initial_setup = true;
579  addExtraVectors();
580 
581  // Perform output related setups
583 
584  // Flush all output to _console that occur during construction and initialization of objects
586 
588  {
589  _resurrector->setRestartFile(_app.getRecoverFileBase());
590  if (_app.getRecoverFileSuffix() == "cpa")
591  _resurrector->setRestartSuffix("xda");
592  }
593 
595  _resurrector->restartFromFile();
596  else
597  {
598  ExodusII_IO * reader = _mesh.exReader();
599 
600  if (reader)
601  {
602  _nl->copyVars(*reader);
603  _aux->copyVars(*reader);
604  }
605  }
606 
607  // Build Refinement and Coarsening maps for stateful material projections if necessary
608  if (_adaptivity.isOn() &&
611  {
613  mooseError("Stateful neighbor material properties do not work with mesh adaptivity");
614 
616  }
617 
618  if (!_app.isRecovering())
619  {
626  {
627  if (!_app.isUltimateMaster())
628  mooseError(
629  "Doing extra refinements when restarting is NOT supported for sub-apps of a MultiApp");
630 
632  }
633  }
634 
635  // Do this just in case things have been done to the mesh
637  _mesh.meshChanged();
638  if (_displaced_problem)
640 
641  unsigned int n_threads = libMesh::n_threads();
642 
643  // UserObject initialSetup
644  std::set<std::string> depend_objects_ic = _ics.getDependObjects();
645  std::set<std::string> depend_objects_aux = _aux->getDependObjects();
646 
647  // This replaces all prior updateDependObjects calls on the old user object warehouses.
648  std::vector<UserObject *> userobjs;
649  theWarehouse().query().condition<AttribSystem>("UserObject").queryInto(userobjs);
650  groupUserObjects(theWarehouse(), userobjs, depend_objects_ic, depend_objects_aux);
651 
652  for (auto obj : userobjs)
653  obj->initialSetup();
654 
655  // check if jacobian calculation is done in userobject
656  for (THREAD_ID tid = 0; tid < n_threads; ++tid)
658  // Check whether nonlocal couling is required or not
662 
663  // Call the initialSetup methods for functions
664  for (THREAD_ID tid = 0; tid < n_threads; tid++)
665  {
667  tid); // initialize scalars so they are properly sized for use as input into ParsedFunctions
669  }
670 
671  // Random interface objects
672  for (const auto & it : _random_data_objects)
673  it.second->updateSeeds(EXEC_INITIAL);
674 
675  if (!_app.isRecovering())
676  {
678 
679  for (THREAD_ID tid = 0; tid < n_threads; tid++)
680  _ics.initialSetup(tid);
682  projectSolution();
683  }
684 
685  // Materials
687  {
688  for (THREAD_ID tid = 0; tid < n_threads; tid++)
689  {
690  // Sort the Material objects, these will be actually computed by MOOSE in reinit methods.
693 
694  // Call initialSetup on both Material and Material objects
696  }
697 
698  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
706  _assembly);
714  cmt(elem_range, true);
715 
719  }
720 
721  for (THREAD_ID tid = 0; tid < n_threads; tid++)
722  {
725  _markers.sort(tid);
726  _markers.initialSetup(tid);
727  }
728 
729 #ifdef LIBMESH_ENABLE_AMR
730 
731  if (!_app.isRecovering())
732  {
733  unsigned int n = adaptivity().getInitialSteps();
734  if (n && !_app.isUltimateMaster() && _app.isRestarting())
735  mooseError("Cannot perform initial adaptivity during restart on sub-apps of a MultiApp!");
736 
738  }
739 
740 #endif // LIBMESH_ENABLE_AMR
741 
742  if (!_app.isRecovering() && !_app.isRestarting())
743  {
744  // During initial setup the solution is copied to solution_old and solution_older
746  }
747 
748  if (!_app.isRecovering())
749  {
750  if (haveXFEM() && updateMeshXFEM())
751  _console << "XFEM updated mesh on initializaton" << std::endl;
752  }
753 
754  // Call initialSetup on the nonlinear system
755  _nl->initialSetup();
756 
757  // Auxilary variable initialSetup calls
758  _aux->initialSetup();
759 
760  _nl->setSolution(*(_nl->system().current_local_solution.get()));
761 
762  // Update the nearest node searches (has to be called after the problem is all set up)
763  // We do this here because this sets up the Element's DoFs to ghost
765 
767  if (_displaced_mesh)
769 
770  // We need to move the mesh in order to build a map between mortar slave and master
771  // interfaces. This map will then be used by the AgumentSparsityOnInterface ghosting functor to
772  // know which dofs we need ghosted when we call EquationSystems::reinit
774  _displaced_problem->updateMesh();
775 
776  // Build the mortar segment meshes for a couple reasons:
777  // 1) Get the ghosting correct for both static and dynamic meshes
778  // 2) Make sure the mortar mesh is built for mortar constraints that live on the static mesh
779  //
780  // It is worth-while to note that mortar meshes that live on a dynamic mesh will be built
781  // during residual and Jacobian evaluation because when displacements are solution variables
782  // the mortar mesh will move and change during the course of a non-linear solve. We DO NOT
783  // redo ghosting during non-linear solve, so for purpose 1) the below call has to be made
785 
786  // Possibly reinit one more time to get ghosting correct
788 
789  if (_displaced_mesh)
790  _displaced_problem->updateMesh();
791 
792  updateGeomSearch(); // Call all of the rest of the geometric searches
793 
794  auto ti = _nl->getTimeIntegrator();
795 
796  if (ti)
797  ti->initialSetup();
798 
799  if (_app.isRestarting() || _app.isRecovering())
800  {
801  if (_app.hasCachedBackup()) // This happens when this app is a sub-app and has been given a
802  // Backup
804  else
805  _resurrector->restartRestartableData();
806 
807  // We may have just clobbered initial conditions that were explicitly set
808  // In a _restart_ scenario it is completely valid to specify new initial conditions
809  // for some of the variables which should override what's coming from the restart file
810  if (!_app.isRecovering())
811  {
812  for (THREAD_ID tid = 0; tid < n_threads; tid++)
813  _ics.initialSetup(tid);
814  _scalar_ics.sort();
815  projectSolution();
816  }
817  }
818 
819  // HUGE NOTE: MultiApp initialSetup() MUST... I repeat MUST be _after_ restartable data has been
820  // restored
821 
822  // Call initialSetup on the MultiApps
823  if (_multi_apps.hasObjects())
824  {
825  _console << COLOR_CYAN << "Initializing All MultiApps" << COLOR_DEFAULT << std::endl;
827  _console << COLOR_CYAN << "Finished Initializing All MultiApps" << COLOR_DEFAULT << std::endl;
828  }
829 
830  // Call initialSetup on the transfers
834 
835  if (!_app.isRecovering())
836  {
838 
840  if (!converged)
841  mooseError("failed to converge initial MultiApp");
842 
843  // We'll backup the Multiapp here
845 
846  for (THREAD_ID tid = 0; tid < n_threads; tid++)
847  reinitScalars(tid);
848 
850 
851  // The FEProblemBase::execute method doesn't call all the systems on EXEC_INITIAL, but it does
852  // set/unset the current flag. Therefore, this resets the current flag to EXEC_INITIAL so that
853  // subsequent calls (e.g., executeControls) have the proper flag.
855  }
856 
857  // Here we will initialize the stateful properties once more since they may have been updated
858  // during initialSetup by calls to computeProperties.
859  //
860  // It's really bad that we don't allow this during restart. It means that we can't add new
861  // stateful materials
862  // during restart. This is only happening because this _has_ to be below initial userobject
863  // execution.
864  // Otherwise this could be done up above... _before_ restoring restartable data... which would
865  // allow you to have
866  // this happen during restart. I honestly have no idea why this has to happen after initial user
867  // object computation.
868  // THAT is something we should fix... so I've opened this ticket: #5804
869  if (!_app.isRecovering() && !_app.isRestarting() &&
872  {
873  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
881  _assembly);
882  Threads::parallel_reduce(elem_range, cmt);
883  }
884 
885  // Control Logic
887 
888  // Scalar variables need to reinited for the initial conditions to be available for output
889  for (unsigned int tid = 0; tid < n_threads; tid++)
890  reinitScalars(tid);
891 
892  if (_displaced_mesh)
893  _displaced_problem->syncSolutions();
894 
895  // Writes all calls to _console from initialSetup() methods
897 
899  {
901  for (THREAD_ID tid = 0; tid < n_threads; ++tid)
902  _assembly[tid]->initNonlocalCoupling();
903  }
904 
907 }
908 
909 void
911 {
912  // Random interface objects
913  for (const auto & it : _random_data_objects)
914  it.second->updateSeeds(EXEC_TIMESTEP_BEGIN);
915 
916  unsigned int n_threads = libMesh::n_threads();
917  for (THREAD_ID tid = 0; tid < n_threads; tid++)
918  {
921  }
922 
923  _aux->timestepSetup();
924  _nl->timestepSetup();
925 
926  for (THREAD_ID tid = 0; tid < n_threads; tid++)
927  {
930  _markers.timestepSetup(tid);
931  }
932 
933  std::vector<UserObject *> userobjs;
934  theWarehouse().query().condition<AttribSystem>("UserObject").queryInto(userobjs);
935  for (auto obj : userobjs)
936  obj->timestepSetup();
937 
938  // Timestep setup of output objects
940 
943  _has_nonlocal_coupling = true;
944 }
945 
946 unsigned int
948 {
949  if (_max_qps == std::numeric_limits<unsigned int>::max())
950  mooseError("Max QPS uninitialized");
951  return _max_qps;
952 }
953 
954 unsigned int
956 {
957  if (_max_shape_funcs == std::numeric_limits<unsigned int>::max())
958  mooseError("Max shape functions uninitialized");
959  return _max_shape_funcs;
960 }
961 
962 Order
964 {
965  return _max_scalar_order;
966 }
967 
968 void
970 {
971  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
972  {
973  const auto & all_kernels = _nl->getKernelWarehouse();
974  const auto & kernels = all_kernels.getObjects(tid);
975  for (const auto & kernel : kernels)
976  {
977  std::shared_ptr<NonlocalKernel> nonlocal_kernel =
979  if (nonlocal_kernel)
980  {
983  _nonlocal_kernels.addObject(kernel, tid);
984  }
985  }
986  const MooseObjectWarehouse<IntegratedBCBase> & all_integrated_bcs =
987  _nl->getIntegratedBCWarehouse();
988  const auto & integrated_bcs = all_integrated_bcs.getObjects(tid);
989  for (const auto & integrated_bc : integrated_bcs)
990  {
991  std::shared_ptr<NonlocalIntegratedBC> nonlocal_integrated_bc =
993  if (nonlocal_integrated_bc)
994  {
997  _nonlocal_integrated_bcs.addObject(integrated_bc, tid);
998  }
999  }
1000  }
1001 }
1002 
1003 void
1005 {
1006  std::set<MooseVariableFEBase *> uo_jacobian_moose_vars;
1007  {
1008  std::vector<ShapeElementUserObject *> objs;
1009  theWarehouse()
1010  .query()
1012  .condition<AttribThread>(tid)
1013  .queryInto(objs);
1014 
1015  for (const auto & uo : objs)
1016  {
1017  _calculate_jacobian_in_uo = uo->computeJacobianFlag();
1018  const std::set<MooseVariableFEBase *> & mv_deps = uo->jacobianMooseVariables();
1019  uo_jacobian_moose_vars.insert(mv_deps.begin(), mv_deps.end());
1020  }
1021  }
1022  {
1023  std::vector<ShapeSideUserObject *> objs;
1024  theWarehouse()
1025  .query()
1027  .condition<AttribThread>(tid)
1028  .queryInto(objs);
1029  for (const auto & uo : objs)
1030  {
1031  _calculate_jacobian_in_uo = uo->computeJacobianFlag();
1032  const std::set<MooseVariableFEBase *> & mv_deps = uo->jacobianMooseVariables();
1033  uo_jacobian_moose_vars.insert(mv_deps.begin(), mv_deps.end());
1034  }
1035  }
1036 
1037  _uo_jacobian_moose_vars[tid].assign(uo_jacobian_moose_vars.begin(), uo_jacobian_moose_vars.end());
1038  std::sort(
1039  _uo_jacobian_moose_vars[tid].begin(), _uo_jacobian_moose_vars[tid].end(), sortMooseVariables);
1040 }
1041 
1042 void
1043 FEProblemBase::setVariableAllDoFMap(const std::vector<MooseVariableFEBase *> moose_vars)
1044 {
1045  for (unsigned int i = 0; i < moose_vars.size(); ++i)
1046  {
1047  VariableName var_name = moose_vars[i]->name();
1048  _nl->setVariableGlobalDoFs(var_name);
1049  _var_dof_map[var_name] = _nl->getVariableGlobalDoFs();
1050  }
1051 }
1052 
1053 void
1054 FEProblemBase::prepare(const Elem * elem, THREAD_ID tid)
1055 {
1056  _assembly[tid]->reinit(elem);
1057 
1058  _nl->prepare(tid);
1059  _aux->prepare(tid);
1060  if (!_has_jacobian || !_const_jacobian)
1061  _assembly[tid]->prepareJacobianBlock();
1062  _assembly[tid]->prepareResidual();
1064  _assembly[tid]->prepareNonlocal();
1065 
1067  {
1068  _displaced_problem->prepare(_displaced_mesh->elemPtr(elem->id()), tid);
1070  _displaced_problem->prepareNonlocal(tid);
1071  }
1072 }
1073 
1074 void
1075 FEProblemBase::prepareFace(const Elem * elem, THREAD_ID tid)
1076 {
1077  _nl->prepareFace(tid, true);
1078  _aux->prepareFace(tid, false);
1079 
1081  _displaced_problem->prepareFace(_displaced_mesh->elemPtr(elem->id()), tid);
1082 }
1083 
1084 void
1085 FEProblemBase::prepare(const Elem * elem,
1086  unsigned int ivar,
1087  unsigned int jvar,
1088  const std::vector<dof_id_type> & dof_indices,
1089  THREAD_ID tid)
1090 {
1091  _assembly[tid]->reinit(elem);
1092 
1093  _nl->prepare(tid);
1094  _aux->prepare(tid);
1095  _assembly[tid]->prepareBlock(ivar, jvar, dof_indices);
1097  if (_nonlocal_cm(ivar, jvar) != 0)
1098  {
1099  MooseVariableFEBase & jv = _nl->getVariable(tid, jvar);
1100  _assembly[tid]->prepareBlockNonlocal(ivar, jvar, dof_indices, jv.allDofIndices());
1101  }
1102 
1104  {
1105  _displaced_problem->prepare(_displaced_mesh->elemPtr(elem->id()), ivar, jvar, dof_indices, tid);
1107  if (_nonlocal_cm(ivar, jvar) != 0)
1108  {
1109  MooseVariableFEBase & jv = _nl->getVariable(tid, jvar);
1110  _displaced_problem->prepareBlockNonlocal(ivar, jvar, dof_indices, jv.allDofIndices(), tid);
1111  }
1112  }
1113 }
1114 
1115 void
1117 {
1118  SubdomainID did = elem->subdomain_id();
1119  _assembly[tid]->setCurrentSubdomainID(did);
1121  _displaced_problem->assembly(tid).setCurrentSubdomainID(did);
1122 }
1123 
1124 void
1125 FEProblemBase::setNeighborSubdomainID(const Elem * elem, unsigned int side, THREAD_ID tid)
1126 {
1127  SubdomainID did = elem->neighbor_ptr(side)->subdomain_id();
1128  _assembly[tid]->setCurrentNeighborSubdomainID(did);
1130  _displaced_problem->assembly(tid).setCurrentNeighborSubdomainID(did);
1131 }
1132 
1133 void
1135 {
1136  SubdomainID did = elem->subdomain_id();
1137  _assembly[tid]->setCurrentNeighborSubdomainID(did);
1139  _displaced_problem->assembly(tid).setCurrentNeighborSubdomainID(did);
1140 }
1141 
1142 void
1144 {
1145  _assembly[tid]->prepare();
1147  _assembly[tid]->prepareNonlocal();
1148 
1150  {
1151  _displaced_problem->prepareAssembly(tid);
1153  _displaced_problem->prepareNonlocal(tid);
1154  }
1155 }
1156 
1157 void
1159 {
1160  _assembly[tid]->addResidual(getVectorTags());
1161 
1162  if (_displaced_problem)
1163  _displaced_problem->addResidual(tid);
1164 }
1165 
1166 void
1168 {
1169  _assembly[tid]->addResidualNeighbor(getVectorTags());
1170 
1171  if (_displaced_problem)
1172  _displaced_problem->addResidualNeighbor(tid);
1173 }
1174 
1175 void
1177 {
1178  _assembly[tid]->addResidualScalar(getVectorTags());
1179 }
1180 
1181 void
1183 {
1184  _assembly[tid]->cacheResidual();
1185  if (_displaced_problem)
1186  _displaced_problem->cacheResidual(tid);
1187 }
1188 
1189 void
1191 {
1192  _assembly[tid]->cacheResidualNeighbor();
1193  if (_displaced_problem)
1194  _displaced_problem->cacheResidualNeighbor(tid);
1195 }
1196 
1197 void
1199 {
1200  _assembly[tid]->addCachedResiduals();
1201 
1202  if (_displaced_problem)
1203  _displaced_problem->addCachedResidual(tid);
1204 }
1205 
1206 void
1207 FEProblemBase::addCachedResidualDirectly(NumericVector<Number> & residual, THREAD_ID tid)
1208 {
1209  _assembly[tid]->addCachedResidual(residual, _nl->timeVectorTag());
1210  _assembly[tid]->addCachedResidual(residual, _nl->nonTimeVectorTag());
1211 
1212  if (_displaced_problem)
1213  _displaced_problem->addCachedResidualDirectly(residual, tid);
1214 }
1215 
1216 void
1217 FEProblemBase::setResidual(NumericVector<Number> & residual, THREAD_ID tid)
1218 {
1219  _assembly[tid]->setResidual(residual);
1220  if (_displaced_problem)
1221  _displaced_problem->setResidual(residual, tid);
1222 }
1223 
1224 void
1225 FEProblemBase::setResidualNeighbor(NumericVector<Number> & residual, THREAD_ID tid)
1226 {
1227  _assembly[tid]->setResidualNeighbor(residual);
1228  if (_displaced_problem)
1229  _displaced_problem->setResidualNeighbor(residual, tid);
1230 }
1231 
1232 void
1234 {
1235  _assembly[tid]->addJacobian();
1237  _assembly[tid]->addJacobianNonlocal();
1238  if (_displaced_problem)
1239  {
1240  _displaced_problem->addJacobian(tid);
1242  _displaced_problem->addJacobianNonlocal(tid);
1243  }
1244 }
1245 
1246 void
1248 {
1249  _assembly[tid]->addJacobianNeighbor();
1250  if (_displaced_problem)
1251  _displaced_problem->addJacobianNeighbor(tid);
1252 }
1253 
1254 void
1256 {
1257  _assembly[tid]->addJacobianScalar();
1258 }
1259 
1260 void
1261 FEProblemBase::addJacobianOffDiagScalar(unsigned int ivar, THREAD_ID tid /* = 0*/)
1262 {
1263  _assembly[tid]->addJacobianOffDiagScalar(ivar);
1264 }
1265 
1266 void
1268 {
1269  _assembly[tid]->cacheJacobian();
1271  _assembly[tid]->cacheJacobianNonlocal();
1272  if (_displaced_problem)
1273  {
1274  _displaced_problem->cacheJacobian(tid);
1276  _displaced_problem->cacheJacobianNonlocal(tid);
1277  }
1278 }
1279 
1280 void
1282 {
1283  _assembly[tid]->cacheJacobianNeighbor();
1284  if (_displaced_problem)
1285  _displaced_problem->cacheJacobianNeighbor(tid);
1286 }
1287 
1288 void
1290 {
1291  _assembly[tid]->addCachedJacobian();
1292  if (_displaced_problem)
1293  _displaced_problem->addCachedJacobian(tid);
1294 }
1295 
1296 void
1297 FEProblemBase::addJacobianBlock(SparseMatrix<Number> & jacobian,
1298  unsigned int ivar,
1299  unsigned int jvar,
1300  const DofMap & dof_map,
1301  std::vector<dof_id_type> & dof_indices,
1302  THREAD_ID tid)
1303 {
1304  _assembly[tid]->addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices);
1306  if (_nonlocal_cm(ivar, jvar) != 0)
1307  {
1308  MooseVariableFEBase & jv = _nl->getVariable(tid, jvar);
1309  _assembly[tid]->addJacobianBlockNonlocal(
1310  jacobian, ivar, jvar, dof_map, dof_indices, jv.allDofIndices());
1311  }
1312 
1313  if (_displaced_problem)
1314  {
1315  _displaced_problem->addJacobianBlock(jacobian, ivar, jvar, dof_map, dof_indices, tid);
1317  if (_nonlocal_cm(ivar, jvar) != 0)
1318  {
1319  MooseVariableFEBase & jv = _nl->getVariable(tid, jvar);
1320  _displaced_problem->addJacobianBlockNonlocal(
1321  jacobian, ivar, jvar, dof_map, dof_indices, jv.allDofIndices(), tid);
1322  }
1323  }
1324 }
1325 
1326 void
1327 FEProblemBase::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
1328  unsigned int ivar,
1329  unsigned int jvar,
1330  const DofMap & dof_map,
1331  std::vector<dof_id_type> & dof_indices,
1332  std::vector<dof_id_type> & neighbor_dof_indices,
1333  THREAD_ID tid)
1334 {
1335  _assembly[tid]->addJacobianNeighbor(
1336  jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices);
1337  if (_displaced_problem)
1338  _displaced_problem->addJacobianNeighbor(
1339  jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices, tid);
1340 }
1341 
1342 void
1344 {
1345  _assembly[tid]->copyShapes(var);
1346 }
1347 
1348 void
1350 {
1351  _assembly[tid]->copyFaceShapes(var);
1352 }
1353 
1354 void
1356 {
1357  _assembly[tid]->copyNeighborShapes(var);
1358 }
1359 
1360 void
1361 FEProblemBase::addGhostedElem(dof_id_type elem_id)
1362 {
1363  if (_mesh.elemPtr(elem_id)->processor_id() != processor_id())
1364  _ghosted_elems.insert(elem_id);
1365 }
1366 
1367 void
1369 {
1370  _mesh.addGhostedBoundary(boundary_id);
1371  if (_displaced_problem)
1372  _displaced_mesh->addGhostedBoundary(boundary_id);
1373 }
1374 
1375 void
1377 {
1379 
1380  if (_displaced_problem)
1382 }
1383 
1384 void
1385 FEProblemBase::sizeZeroes(unsigned int /*size*/, THREAD_ID /*tid*/)
1386 {
1387  mooseDoOnce(mooseWarning(
1388  "This function is deprecated and no longer performs any function. Please do not call it."));
1389 }
1390 
1391 bool
1392 FEProblemBase::reinitDirac(const Elem * elem, THREAD_ID tid)
1393 {
1394  std::vector<Point> & points = _dirac_kernel_info.getPoints()[elem].first;
1395 
1396  unsigned int n_points = points.size();
1397 
1398  if (n_points)
1399  {
1400  if (n_points > _max_qps)
1401  {
1402  _max_qps = n_points;
1403 
1408  unsigned int max_qpts = getMaxQps();
1409  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
1410  {
1411  // the highest available order in libMesh is 43
1412  _scalar_zero[tid].resize(FORTYTHIRD, 0);
1413  _zero[tid].resize(max_qpts, 0);
1414  _grad_zero[tid].resize(max_qpts, RealGradient(0.));
1415  _second_zero[tid].resize(max_qpts, RealTensor(0.));
1416  _second_phi_zero[tid].resize(
1417  max_qpts, std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));
1418  _vector_zero[tid].resize(max_qpts, RealGradient(0.));
1419  _vector_curl_zero[tid].resize(max_qpts, RealGradient(0.));
1420  }
1421  }
1422 
1423  _assembly[tid]->reinitAtPhysical(elem, points);
1424 
1425  _nl->prepare(tid);
1426  _aux->prepare(tid);
1427 
1428  reinitElem(elem, tid);
1429  }
1430 
1431  _assembly[tid]->prepare();
1433  _assembly[tid]->prepareNonlocal();
1434 
1435  bool have_points = n_points > 0;
1437  {
1438  have_points |= _displaced_problem->reinitDirac(_displaced_mesh->elemPtr(elem->id()), tid);
1440  _displaced_problem->prepareNonlocal(tid);
1441  }
1442 
1443  return have_points;
1444 }
1445 
1446 void
1447 FEProblemBase::reinitElem(const Elem * elem, THREAD_ID tid)
1448 {
1449  _nl->reinitElem(elem, tid);
1450  _aux->reinitElem(elem, tid);
1451 
1453  _displaced_problem->reinitElem(_displaced_mesh->elemPtr(elem->id()), tid);
1454 }
1455 
1456 void
1458  const std::vector<Point> & phys_points_in_elem,
1459  THREAD_ID tid,
1460  bool suppress_displaced_init)
1461 {
1462  _assembly[tid]->reinitAtPhysical(elem, phys_points_in_elem);
1463 
1464  _nl->prepare(tid);
1465  _aux->prepare(tid);
1466 
1467  reinitElem(elem, tid);
1468  _assembly[tid]->prepare();
1470  _assembly[tid]->prepareNonlocal();
1471 
1472  if (_displaced_problem && _reinit_displaced_elem && !suppress_displaced_init)
1473  {
1474  _displaced_problem->reinitElemPhys(
1475  _displaced_mesh->elemPtr(elem->id()), phys_points_in_elem, tid);
1477  _displaced_problem->prepareNonlocal(tid);
1478  }
1479 }
1480 
1481 void
1483  unsigned int side,
1484  BoundaryID bnd_id,
1485  THREAD_ID tid)
1486 {
1487  _assembly[tid]->reinit(elem, side);
1488 
1489  _nl->reinitElemFace(elem, side, bnd_id, tid);
1490  _aux->reinitElemFace(elem, side, bnd_id, tid);
1491 
1493  _displaced_problem->reinitElemFace(_displaced_mesh->elemPtr(elem->id()), side, bnd_id, tid);
1494 }
1495 
1496 void
1497 FEProblemBase::reinitNode(const Node * node, THREAD_ID tid)
1498 {
1499  _assembly[tid]->reinit(node);
1500 
1502  _displaced_problem->reinitNode(&_displaced_mesh->nodeRef(node->id()), tid);
1503 
1504  _nl->reinitNode(node, tid);
1505  _aux->reinitNode(node, tid);
1506 }
1507 
1508 void
1509 FEProblemBase::reinitNodeFace(const Node * node, BoundaryID bnd_id, THREAD_ID tid)
1510 {
1511  _assembly[tid]->reinit(node);
1512 
1514  _displaced_problem->reinitNodeFace(&_displaced_mesh->nodeRef(node->id()), bnd_id, tid);
1515 
1516  _nl->reinitNodeFace(node, bnd_id, tid);
1517  _aux->reinitNodeFace(node, bnd_id, tid);
1518 }
1519 
1520 void
1521 FEProblemBase::reinitNodes(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
1522 {
1524  _displaced_problem->reinitNodes(nodes, tid);
1525 
1526  _nl->reinitNodes(nodes, tid);
1527  _aux->reinitNodes(nodes, tid);
1528 }
1529 
1530 void
1531 FEProblemBase::reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, THREAD_ID tid)
1532 {
1534  _displaced_problem->reinitNodesNeighbor(nodes, tid);
1535 
1536  _nl->reinitNodesNeighbor(nodes, tid);
1537  _aux->reinitNodesNeighbor(nodes, tid);
1538 }
1539 
1540 void
1542 {
1544  _displaced_problem->reinitScalars(tid);
1545 
1546  _nl->reinitScalars(tid);
1547  _aux->reinitScalars(tid);
1548 
1549  _assembly[tid]->prepareScalar();
1550 }
1551 
1552 void
1554 {
1555  _assembly[tid]->prepareOffDiagScalar();
1556  if (_displaced_problem)
1557  _displaced_problem->reinitOffDiagScalars(tid);
1558 }
1559 
1560 void
1561 FEProblemBase::reinitNeighbor(const Elem * elem, unsigned int side, THREAD_ID tid)
1562 {
1563  const Elem * neighbor = elem->neighbor_ptr(side);
1564  unsigned int neighbor_side = neighbor->which_neighbor_am_i(elem);
1565 
1566  _assembly[tid]->reinitElemAndNeighbor(elem, side, neighbor, neighbor_side);
1567 
1568  _nl->prepareNeighbor(tid);
1569  _aux->prepareNeighbor(tid);
1570 
1571  _assembly[tid]->prepareNeighbor();
1572 
1573  BoundaryID bnd_id = 0; // some dummy number (it is not really used for anything, right now)
1574  _nl->reinitElemFace(elem, side, bnd_id, tid);
1575  _aux->reinitElemFace(elem, side, bnd_id, tid);
1576 
1577  _nl->reinitNeighborFace(neighbor, neighbor_side, bnd_id, tid);
1578  _aux->reinitNeighborFace(neighbor, neighbor_side, bnd_id, tid);
1579 
1581  _displaced_problem->reinitNeighbor(elem, side, tid);
1582 }
1583 
1584 void
1586  unsigned int neighbor_side,
1587  const std::vector<Point> & physical_points,
1588  THREAD_ID tid)
1589 {
1590  // Reinits shape the functions at the physical points
1591  _assembly[tid]->reinitNeighborAtPhysical(neighbor, neighbor_side, physical_points);
1592 
1593  // Sets the neighbor dof indices
1594  _nl->prepareNeighbor(tid);
1595  _aux->prepareNeighbor(tid);
1596 
1597  // Resizes Re and Ke
1598  _assembly[tid]->prepareNeighbor();
1599 
1600  // Compute the values of each variable at the points
1601  _nl->reinitNeighborFace(neighbor, neighbor_side, 0, tid);
1602  _aux->reinitNeighborFace(neighbor, neighbor_side, 0, tid);
1603 
1604  // Do the same for the displaced problem
1606  _displaced_problem->reinitNeighborPhys(
1607  _displaced_mesh->elemPtr(neighbor->id()), neighbor_side, physical_points, tid);
1608 }
1609 
1610 void
1612  const std::vector<Point> & physical_points,
1613  THREAD_ID tid)
1614 {
1615  // Reinits shape the functions at the physical points
1616  _assembly[tid]->reinitNeighborAtPhysical(neighbor, physical_points);
1617 
1618  // Sets the neighbor dof indices
1619  _nl->prepareNeighbor(tid);
1620  _aux->prepareNeighbor(tid);
1621 
1622  // Resizes Re and Ke
1623  _assembly[tid]->prepareNeighbor();
1624 
1625  // Compute the values of each variable at the points
1626  _nl->reinitNeighbor(neighbor, tid);
1627  _aux->reinitNeighbor(neighbor, tid);
1628 
1629  // Do the same for the displaced problem
1631  _displaced_problem->reinitNeighborPhys(
1632  _displaced_mesh->elemPtr(neighbor->id()), physical_points, tid);
1633 }
1634 
1635 void
1636 FEProblemBase::getDiracElements(std::set<const Elem *> & elems)
1637 {
1638  // First add in the undisplaced elements
1639  elems = _dirac_kernel_info.getElements();
1640 
1641  if (_displaced_problem)
1642  {
1643  std::set<const Elem *> displaced_elements;
1644  _displaced_problem->getDiracElements(displaced_elements);
1645 
1646  { // Use the ids from the displaced elements to get the undisplaced elements
1647  // and add them to the list
1648  for (const auto & elem : displaced_elements)
1649  elems.insert(_mesh.elemPtr(elem->id()));
1650  }
1651  }
1652 }
1653 
1654 void
1656 {
1658 
1659  if (_displaced_problem)
1660  _displaced_problem->clearDiracInfo();
1661 }
1662 
1663 void
1665 {
1666  _all_materials.subdomainSetup(subdomain, tid);
1667 
1668  // Call the subdomain methods of the output system, these are not threaded so only call it once
1669  if (tid == 0)
1671 
1672  _nl->subdomainSetup(subdomain, tid);
1673 
1674  // FIXME: call displaced_problem->subdomainSetup() ?
1675  // When adding possibility with materials being evaluated on displaced mesh
1676 }
1677 
1678 void
1680 {
1681  _all_materials.neighborSubdomainSetup(subdomain, tid);
1682 }
1683 
1684 void
1685 FEProblemBase::addFunction(std::string type, const std::string & name, InputParameters parameters)
1686 {
1687  parameters.set<SubProblem *>("_subproblem") = this;
1688 
1689  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1690  {
1691  std::shared_ptr<Function> func = _factory.create<Function>(type, name, parameters, tid);
1692  _functions.addObject(func, tid);
1693  }
1694 }
1695 
1696 bool
1697 FEProblemBase::hasFunction(const std::string & name, THREAD_ID tid)
1698 {
1699  return _functions.hasActiveObject(name, tid);
1700 }
1701 
1702 Function &
1703 FEProblemBase::getFunction(const std::string & name, THREAD_ID tid)
1704 {
1705  // This thread lock is necessary since this method will create functions
1706  // for all threads if one is missing.
1707  Threads::spin_mutex::scoped_lock lock(get_function_mutex);
1708 
1709  if (!hasFunction(name, tid))
1710  {
1711  // If we didn't find a function, it might be a default function, attempt to construct one now
1712  std::istringstream ss(name);
1713  Real real_value;
1714 
1715  // First see if it's just a constant. If it is, build a ConstantFunction
1716  if (ss >> real_value && ss.eof())
1717  {
1718  InputParameters params = _factory.getValidParams("ConstantFunction");
1719  params.set<Real>("value") = real_value;
1720  addFunction("ConstantFunction", ss.str(), params);
1721  }
1722  else
1723  {
1725  std::string vars = "x,y,z,t,NaN,pi,e";
1726  if (fp.Parse(name, vars) == -1) // -1 for success
1727  {
1728  // It parsed ok, so build a MooseParsedFunction
1729  InputParameters params = _factory.getValidParams("ParsedFunction");
1730  params.set<std::string>("value") = name;
1731  addFunction("ParsedFunction", name, params);
1732  }
1733  }
1734 
1735  // Try once more
1736  if (!hasFunction(name, tid))
1737  mooseError("Unable to find function " + name);
1738  }
1739 
1740  return *(_functions.getActiveObject(name, tid));
1741 }
1742 
1743 void
1745 {
1746  _line_search->lineSearch();
1747 }
1748 
1751 {
1752  mooseDeprecated("FEProblemBase::getNonlinearSystem() is deprecated, please use "
1753  "FEProblemBase::getNonlinearSystemBase() \n");
1754 
1756 
1757  if (!nl_sys)
1758  mooseError("This is not a NonlinearSystem");
1759 
1760  return *nl_sys;
1761 }
1762 
1763 void
1765  const std::string & name,
1766  InputParameters parameters)
1767 {
1768  parameters.set<std::string>("type") = type;
1769  std::shared_ptr<Distribution> dist = _factory.create<Distribution>(type, name, parameters);
1770  _distributions.addObject(dist);
1771 }
1772 
1773 Distribution &
1774 FEProblemBase::getDistribution(const std::string & name)
1775 {
1777  mooseError("Unable to find distribution " + name);
1778 
1779  return *(_distributions.getActiveObject(name));
1780 }
1781 
1782 void
1783 FEProblemBase::addSampler(std::string type, const std::string & name, InputParameters parameters)
1784 {
1785  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
1786  {
1787  std::shared_ptr<Sampler> dist = _factory.create<Sampler>(type, name, parameters, tid);
1788  _samplers.addObject(dist, tid);
1789  }
1790 }
1791 
1792 Sampler &
1793 FEProblemBase::getSampler(const std::string & name, THREAD_ID tid)
1794 {
1795  if (!_samplers.hasActiveObject(name, tid))
1796  mooseError("Unable to find Sampler " + name);
1797 
1798  return *(_samplers.getActiveObject(name, tid));
1799 }
1800 
1801 bool
1802 FEProblemBase::duplicateVariableCheck(const std::string & var_name,
1803  const FEType & type,
1804  bool is_aux)
1805 {
1806  SystemBase * curr_sys_ptr = _nl.get();
1807  SystemBase * other_sys_ptr = _aux.get();
1808  std::string error_prefix = "";
1809  if (is_aux)
1810  {
1811  curr_sys_ptr = _aux.get();
1812  other_sys_ptr = _nl.get();
1813  error_prefix = "Aux";
1814  }
1815 
1816  if (other_sys_ptr->hasVariable(var_name))
1817  mooseError("Cannot have an auxiliary variable and a nonlinear variable with the same name: ",
1818  var_name);
1819 
1820  if (curr_sys_ptr->hasVariable(var_name))
1821  {
1822  const Variable & var =
1823  curr_sys_ptr->system().variable(curr_sys_ptr->system().variable_number(var_name));
1824  if (var.type() != type)
1825  mooseError(error_prefix,
1826  "Variable with name '",
1827  var_name,
1828  "' already exists but is of a differing type!");
1829 
1830  return true;
1831  }
1832 
1833  return false;
1834 }
1835 
1836 void
1837 FEProblemBase::addVariable(const std::string & var_name,
1838  const FEType & type,
1839  Real scale_factor,
1840  const std::set<SubdomainID> * const active_subdomains)
1841 {
1842  if (duplicateVariableCheck(var_name, type, /* is_aux = */ false))
1843  return;
1844 
1845  _nl->addVariable(var_name, type, scale_factor, active_subdomains);
1846  if (_displaced_problem)
1847  _displaced_problem->addVariable(var_name, type, scale_factor, active_subdomains);
1848 }
1849 
1850 void
1851 FEProblemBase::addScalarVariable(const std::string & var_name,
1852  Order order,
1853  Real scale_factor,
1854  const std::set<SubdomainID> * const active_subdomains)
1855 {
1856  if (order > _max_scalar_order)
1857  _max_scalar_order = order;
1858 
1859  FEType type(order, SCALAR);
1860  if (duplicateVariableCheck(var_name, type, /* is_aux = */ false))
1861  return;
1862 
1863  _nl->addScalarVariable(var_name, order, scale_factor, active_subdomains);
1864  if (_displaced_problem)
1865  _displaced_problem->addScalarVariable(var_name, order, scale_factor, active_subdomains);
1866 }
1867 
1868 void
1869 FEProblemBase::addKernel(const std::string & kernel_name,
1870  const std::string & name,
1871  InputParameters parameters)
1872 {
1873  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
1874  {
1875  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1876  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1877  const auto & disp_names = _displaced_problem->getDisplacementVarNames();
1878  parameters.set<std::vector<VariableName>>("displacements") =
1879  std::vector<VariableName>(disp_names.begin(), disp_names.end());
1880  _reinit_displaced_elem = true;
1881  }
1882  else
1883  {
1884  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
1885  {
1886  // We allow Kernels to request that they use_displaced_mesh,
1887  // but then be overridden when no displacements variables are
1888  // provided in the Mesh block. If that happened, update the value
1889  // of use_displaced_mesh appropriately for this Kernel.
1890  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1891  parameters.set<bool>("use_displaced_mesh") = false;
1892  }
1893 
1894  parameters.set<SubProblem *>("_subproblem") = this;
1895  parameters.set<SystemBase *>("_sys") = _nl.get();
1896  }
1897 
1898  _nl->addKernel(kernel_name, name, parameters);
1899 }
1900 
1901 void
1902 FEProblemBase::addNodalKernel(const std::string & kernel_name,
1903  const std::string & name,
1904  InputParameters parameters)
1905 {
1906  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
1907  {
1908  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1909  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1910  _reinit_displaced_elem = true;
1911  }
1912  else
1913  {
1914  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
1915  {
1916  // We allow NodalKernels to request that they use_displaced_mesh,
1917  // but then be overridden when no displacements variables are
1918  // provided in the Mesh block. If that happened, update the value
1919  // of use_displaced_mesh appropriately for this NodalKernel.
1920  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1921  parameters.set<bool>("use_displaced_mesh") = false;
1922  }
1923 
1924  parameters.set<SubProblem *>("_subproblem") = this;
1925  parameters.set<SystemBase *>("_sys") = _nl.get();
1926  }
1927  _nl->addNodalKernel(kernel_name, name, parameters);
1928 }
1929 
1930 void
1931 FEProblemBase::addScalarKernel(const std::string & kernel_name,
1932  const std::string & name,
1933  InputParameters parameters)
1934 {
1935  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
1936  {
1937  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1938  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1939  }
1940  else
1941  {
1942  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
1943  {
1944  // We allow ScalarKernels to request that they use_displaced_mesh,
1945  // but then be overridden when no displacements variables are
1946  // provided in the Mesh block. If that happened, update the value
1947  // of use_displaced_mesh appropriately for this ScalarKernel.
1948  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1949  parameters.set<bool>("use_displaced_mesh") = false;
1950  }
1951 
1952  parameters.set<SubProblem *>("_subproblem") = this;
1953  parameters.set<SystemBase *>("_sys") = _nl.get();
1954  }
1955 
1956  _nl->addScalarKernel(kernel_name, name, parameters);
1957 }
1958 
1959 void
1960 FEProblemBase::addBoundaryCondition(const std::string & bc_name,
1961  const std::string & name,
1962  InputParameters parameters)
1963 {
1964  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
1965  {
1966  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
1967  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
1968  const auto & disp_names = _displaced_problem->getDisplacementVarNames();
1969  parameters.set<std::vector<VariableName>>("displacements") =
1970  std::vector<VariableName>(disp_names.begin(), disp_names.end());
1971  _reinit_displaced_face = true;
1972  }
1973  else
1974  {
1975  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
1976  {
1977  // We allow Materials to request that they use_displaced_mesh,
1978  // but then be overridden when no displacements variables are
1979  // provided in the Mesh block. If that happened, update the value
1980  // of use_displaced_mesh appropriately for this Material.
1981  if (parameters.have_parameter<bool>("use_displaced_mesh"))
1982  parameters.set<bool>("use_displaced_mesh") = false;
1983  }
1984 
1985  parameters.set<SubProblem *>("_subproblem") = this;
1986  parameters.set<SystemBase *>("_sys") = _nl.get();
1987  }
1988 
1989  _nl->addBoundaryCondition(bc_name, name, parameters);
1990 }
1991 
1992 void
1993 FEProblemBase::addConstraint(const std::string & c_name,
1994  const std::string & name,
1995  InputParameters parameters)
1996 {
1997  _has_constraints = true;
1998 
1999  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2000  {
2001  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2002  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2003  _reinit_displaced_face = true;
2004  }
2005  else
2006  {
2007  // It might _want_ to use a displaced mesh... but we're not so set it to false
2008  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2009  parameters.set<bool>("use_displaced_mesh") = false;
2010 
2011  parameters.set<SubProblem *>("_subproblem") = this;
2012  parameters.set<SystemBase *>("_sys") = _nl.get();
2013  }
2014 
2015  // Check that "variable" is in the NonlinearSystem.
2016  if (!_nl->hasVariable(parameters.get<NonlinearVariableName>("variable")))
2017  mooseError(name,
2018  ": Cannot add Constraint for variable ",
2019  parameters.get<NonlinearVariableName>("variable"),
2020  ", it is not a nonlinear variable!");
2021 
2022  _nl->addConstraint(c_name, name, parameters);
2023 }
2024 
2025 void
2026 FEProblemBase::addAuxVariable(const std::string & var_name,
2027  const FEType & type,
2028  const std::set<SubdomainID> * const active_subdomains)
2029 {
2030  if (duplicateVariableCheck(var_name, type, /* is_aux = */ true))
2031  return;
2032 
2033  _aux->addVariable(var_name, type, 1.0, active_subdomains);
2034  if (_displaced_problem)
2035  _displaced_problem->addAuxVariable(var_name, type, active_subdomains);
2036 }
2037 
2038 void
2039 FEProblemBase::addAuxScalarVariable(const std::string & var_name,
2040  Order order,
2041  Real scale_factor,
2042  const std::set<SubdomainID> * const active_subdomains)
2043 {
2044  if (order > _max_scalar_order)
2045  _max_scalar_order = order;
2046 
2047  FEType type(order, SCALAR);
2048  if (duplicateVariableCheck(var_name, type, /* is_aux = */ true))
2049  return;
2050 
2051  _aux->addScalarVariable(var_name, order, scale_factor, active_subdomains);
2052  if (_displaced_problem)
2053  _displaced_problem->addAuxScalarVariable(var_name, order, scale_factor, active_subdomains);
2054 }
2055 
2056 void
2057 FEProblemBase::addAuxKernel(const std::string & kernel_name,
2058  const std::string & name,
2059  InputParameters parameters)
2060 {
2061  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2062  {
2063  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2064  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
2065  parameters.set<SystemBase *>("_nl_sys") = &_displaced_problem->nlSys();
2066  if (!parameters.get<std::vector<BoundaryName>>("boundary").empty())
2067  _reinit_displaced_face = true;
2068  else
2069  _reinit_displaced_elem = true;
2070  }
2071  else
2072  {
2073  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2074  {
2075  // We allow AuxKernels to request that they use_displaced_mesh,
2076  // but then be overridden when no displacements variables are
2077  // provided in the Mesh block. If that happened, update the value
2078  // of use_displaced_mesh appropriately for this AuxKernel.
2079  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2080  parameters.set<bool>("use_displaced_mesh") = false;
2081  }
2082 
2083  parameters.set<SubProblem *>("_subproblem") = this;
2084  parameters.set<SystemBase *>("_sys") = _aux.get();
2085  parameters.set<SystemBase *>("_nl_sys") = _nl.get();
2086  }
2087 
2088  _aux->addKernel(kernel_name, name, parameters);
2089 }
2090 
2091 void
2092 FEProblemBase::addAuxScalarKernel(const std::string & kernel_name,
2093  const std::string & name,
2094  InputParameters parameters)
2095 {
2096  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2097  {
2098  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2099  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
2100  }
2101  else
2102  {
2103  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2104  {
2105  // We allow AuxScalarKernels to request that they use_displaced_mesh,
2106  // but then be overridden when no displacements variables are
2107  // provided in the Mesh block. If that happened, update the value
2108  // of use_displaced_mesh appropriately for this AuxScalarKernel.
2109  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2110  parameters.set<bool>("use_displaced_mesh") = false;
2111  }
2112 
2113  parameters.set<SubProblem *>("_subproblem") = this;
2114  parameters.set<SystemBase *>("_sys") = _aux.get();
2115  }
2116 
2117  _aux->addScalarKernel(kernel_name, name, parameters);
2118 }
2119 
2120 void
2121 FEProblemBase::addDiracKernel(const std::string & kernel_name,
2122  const std::string & name,
2123  InputParameters parameters)
2124 {
2125  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2126  {
2127  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2128  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2129  _reinit_displaced_elem = true;
2130  }
2131  else
2132  {
2133  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2134  {
2135  // We allow DiracKernels to request that they use_displaced_mesh,
2136  // but then be overridden when no displacements variables are
2137  // provided in the Mesh block. If that happened, update the value
2138  // of use_displaced_mesh appropriately for this DiracKernel.
2139  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2140  parameters.set<bool>("use_displaced_mesh") = false;
2141  }
2142 
2143  parameters.set<SubProblem *>("_subproblem") = this;
2144  parameters.set<SystemBase *>("_sys") = _nl.get();
2145  }
2146 
2147  _nl->addDiracKernel(kernel_name, name, parameters);
2148 }
2149 
2150 // DGKernels ////
2151 
2152 void
2153 FEProblemBase::addDGKernel(const std::string & dg_kernel_name,
2154  const std::string & name,
2155  InputParameters parameters)
2156 {
2157  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2158  {
2159  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2160  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2161  _reinit_displaced_face = true;
2162  }
2163  else
2164  {
2165  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2166  {
2167  // We allow DGKernels to request that they use_displaced_mesh,
2168  // but then be overridden when no displacements variables are
2169  // provided in the Mesh block. If that happened, update the value
2170  // of use_displaced_mesh appropriately for this DGKernel.
2171  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2172  parameters.set<bool>("use_displaced_mesh") = false;
2173  }
2174 
2175  parameters.set<SubProblem *>("_subproblem") = this;
2176  parameters.set<SystemBase *>("_sys") = _nl.get();
2177  }
2178 
2179  _nl->addDGKernel(dg_kernel_name, name, parameters);
2180 
2182 }
2183 
2184 // InterfaceKernels ////
2185 
2186 void
2187 FEProblemBase::addInterfaceKernel(const std::string & interface_kernel_name,
2188  const std::string & name,
2189  InputParameters parameters)
2190 {
2191  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2192  {
2193  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2194  parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys();
2195  _reinit_displaced_face = true;
2196  }
2197  else
2198  {
2199  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2200  {
2201  // We allow InterfaceKernels to request that they use_displaced_mesh,
2202  // but then be overridden when no displacements variables are
2203  // provided in the Mesh block. If that happened, update the value
2204  // of use_displaced_mesh appropriately for this InterfaceKernel.
2205  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2206  parameters.set<bool>("use_displaced_mesh") = false;
2207  }
2208 
2209  parameters.set<SubProblem *>("_subproblem") = this;
2210  parameters.set<SystemBase *>("_sys") = _nl.get();
2211  }
2212 
2213  _nl->addInterfaceKernel(interface_kernel_name, name, parameters);
2214 
2216 }
2217 
2218 void
2219 FEProblemBase::addInitialCondition(const std::string & ic_name,
2220  const std::string & name,
2221  InputParameters parameters)
2222 {
2223 
2224  // before we start to mess with the initial condition, we need to check parameters for errors.
2226 
2227  parameters.set<SubProblem *>("_subproblem") = this;
2228 
2229  const std::string & var_name = parameters.get<VariableName>("variable");
2230  // field IC
2231  if (hasVariable(var_name))
2232  {
2233  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2234  {
2237  parameters.set<SystemBase *>("_sys") = &var.sys();
2238  std::shared_ptr<InitialConditionBase> ic;
2239  if (dynamic_cast<MooseVariable *>(&var))
2240  ic = _factory.create<InitialCondition>(ic_name, name, parameters, tid);
2241  else if (dynamic_cast<VectorMooseVariable *>(&var))
2242  ic = _factory.create<VectorInitialCondition>(ic_name, name, parameters, tid);
2243  else
2244  mooseError("Your FE variable in initial condition ",
2245  name,
2246  " must be either of scalar or vector type");
2247  _ics.addObject(ic, tid);
2248  }
2249  }
2250 
2251  // scalar IC
2252  else if (hasScalarVariable(var_name))
2253  {
2254  MooseVariableScalar & var = getScalarVariable(0, var_name);
2255  parameters.set<SystemBase *>("_sys") = &var.sys();
2256  std::shared_ptr<ScalarInitialCondition> ic =
2258  _scalar_ics.addObject(ic);
2259  }
2260 
2261  else
2262  mooseError(
2263  "Variable '", var_name, "' requested in initial condition '", name, "' does not exist.");
2264 }
2265 
2266 void
2268 {
2269  TIME_SECTION(_project_solution_timer)
2270 
2271  FloatingPointExceptionGuard fpe_guard(_app);
2272 
2273  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
2274  ComputeInitialConditionThread cic(*this);
2275  Threads::parallel_reduce(elem_range, cic);
2276 
2277  // Need to close the solution vector here so that boundary ICs take precendence
2278  _nl->solution().close();
2279  _aux->solution().close();
2280 
2281  // now run boundary-restricted initial conditions
2282  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
2284  Threads::parallel_reduce(bnd_nodes, cbic);
2285 
2286  _nl->solution().close();
2287  _aux->solution().close();
2288 
2289  // Also, load values into the SCALAR dofs
2290  // Note: We assume that all SCALAR dofs are on the
2291  // processor with highest ID
2292  if (processor_id() == (n_processors() - 1) && _scalar_ics.hasActiveObjects())
2293  {
2294  const auto & ics = _scalar_ics.getActiveObjects();
2295  for (const auto & ic : ics)
2296  {
2297  MooseVariableScalar & var = ic->variable();
2298  var.reinit();
2299 
2300  DenseVector<Number> vals(var.order());
2301  ic->compute(vals);
2302 
2303  const unsigned int n_SCALAR_dofs = var.dofIndices().size();
2304  for (unsigned int i = 0; i < n_SCALAR_dofs; i++)
2305  {
2306  const dof_id_type global_index = var.dofIndices()[i];
2307  var.sys().solution().set(global_index, vals(i));
2308  var.setValue(i, vals(i));
2309  }
2310  }
2311  }
2312 
2313  _nl->solution().close();
2314  _nl->solution().localize(*_nl->system().current_local_solution, _nl->dofMap().get_send_list());
2315 
2316  _aux->solution().close();
2317  _aux->solution().localize(*_aux->sys().current_local_solution, _aux->dofMap().get_send_list());
2318 }
2319 
2320 std::shared_ptr<Material>
2323  THREAD_ID tid,
2324  bool no_warn)
2325 {
2326  switch (type)
2327  {
2329  name += "_neighbor";
2330  break;
2332  name += "_face";
2333  break;
2334  default:
2335  break;
2336  }
2337 
2338  std::shared_ptr<Material> material = _all_materials[type].getActiveObject(name, tid);
2339  if (!no_warn && material->getParam<bool>("compute") && type == Moose::BLOCK_MATERIAL_DATA)
2340  mooseWarning("You are retrieving a Material object (",
2341  material->name(),
2342  "), but its compute flag is set to true. This indicates that MOOSE is "
2343  "computing this property which may not be desired and produce un-expected "
2344  "results.");
2345 
2346  return material;
2347 }
2348 
2349 std::shared_ptr<MaterialData>
2351 {
2352  std::shared_ptr<MaterialData> output;
2353  switch (type)
2354  {
2356  output = _material_data[tid];
2357  break;
2359  output = _neighbor_material_data[tid];
2360  break;
2363  output = _bnd_material_data[tid];
2364  break;
2365  }
2366  return output;
2367 }
2368 
2369 void
2370 FEProblemBase::addMaterial(const std::string & mat_name,
2371  const std::string & name,
2372  InputParameters parameters)
2373 {
2375 }
2376 
2377 void
2378 FEProblemBase::addADResidualMaterial(const std::string & mat_name,
2379  const std::string & name,
2380  InputParameters parameters)
2381 {
2383 }
2384 
2385 void
2386 FEProblemBase::addADJacobianMaterial(const std::string & mat_name,
2387  const std::string & name,
2388  InputParameters parameters)
2389 {
2391 }
2392 
2393 void
2394 FEProblemBase::addMaterialHelper(std::vector<MaterialWarehouse *> warehouses,
2395  const std::string & mat_name,
2396  const std::string & name,
2397  InputParameters parameters)
2398 {
2399  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2400  {
2401  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2402  _reinit_displaced_elem = true;
2403  }
2404  else
2405  {
2406  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2407  {
2408  // We allow Materials to request that they use_displaced_mesh,
2409  // but then be overridden when no displacements variables are
2410  // provided in the Mesh block. If that happened, update the value
2411  // of use_displaced_mesh appropriately for this Material.
2412  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2413  parameters.set<bool>("use_displaced_mesh") = false;
2414  }
2415 
2416  parameters.set<SubProblem *>("_subproblem") = this;
2417  }
2418 
2419  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
2420  {
2421  // Create the general Block/Boundary Material object
2422  std::shared_ptr<Material> material = _factory.create<Material>(mat_name, name, parameters, tid);
2423  bool discrete = !material->getParam<bool>("compute");
2424 
2425  // If the object is boundary restricted do not create the neighbor and face objects
2426  if (material->boundaryRestricted())
2427  {
2428  _all_materials.addObject(material, tid);
2429  if (discrete)
2430  _discrete_materials.addObject(material, tid);
2431  else
2432  for (auto && warehouse : warehouses)
2433  warehouse->addObject(material, tid);
2434  }
2435 
2436  // Non-boundary restricted require face and neighbor objects
2437  else
2438  {
2439  // The name of the object being created, this is changed multiple times as objects are created
2440  // below
2441  std::string object_name;
2442 
2443  // Create a copy of the supplied parameters to the setting for "_material_data_type" isn't
2444  // used from a previous tid loop
2445  InputParameters current_parameters = parameters;
2446 
2447  // face material
2448  current_parameters.set<Moose::MaterialDataType>("_material_data_type") =
2450  object_name = name + "_face";
2451  std::shared_ptr<Material> face_material =
2452  _factory.create<Material>(mat_name, object_name, current_parameters, tid);
2453 
2454  // neighbor material
2455  current_parameters.set<Moose::MaterialDataType>("_material_data_type") =
2457  current_parameters.set<bool>("_neighbor") = true;
2458  object_name = name + "_neighbor";
2459  std::shared_ptr<Material> neighbor_material =
2460  _factory.create<Material>(mat_name, object_name, current_parameters, tid);
2461 
2462  // Store the material objects
2463  _all_materials.addObjects(material, neighbor_material, face_material, tid);
2464 
2465  if (discrete)
2466  _discrete_materials.addObjects(material, neighbor_material, face_material, tid);
2467  else
2468  for (auto && warehouse : warehouses)
2469  warehouse->addObjects(material, neighbor_material, face_material, tid);
2470 
2471  // link parameters of face and neighbor materials
2472  MooseObjectParameterName name(MooseObjectName("Material", material->name()), "*");
2473  MooseObjectParameterName face_name(MooseObjectName("Material", face_material->name()), "*");
2474  MooseObjectParameterName neighbor_name(MooseObjectName("Material", neighbor_material->name()),
2475  "*");
2478  name, neighbor_name, false);
2479  }
2480  }
2481 }
2482 
2483 void
2485 {
2486  std::set<MooseVariableFEBase *> needed_moose_vars;
2487  std::set<unsigned int> needed_mat_props;
2488 
2489  if (_all_materials.hasActiveBlockObjects(blk_id, tid))
2490  {
2491  _all_materials.updateVariableDependency(needed_moose_vars, tid);
2492  _all_materials.updateBlockMatPropDependency(blk_id, needed_mat_props, tid);
2493  }
2494 
2495  const std::set<BoundaryID> & ids = _mesh.getSubdomainBoundaryIds(blk_id);
2496  for (const auto & id : ids)
2497  {
2499  {
2500  _jacobian_materials.updateBoundaryVariableDependency(id, needed_moose_vars, tid);
2501  _jacobian_materials.updateBoundaryMatPropDependency(id, needed_mat_props, tid);
2502  }
2503  else
2504  {
2505  _residual_materials.updateBoundaryVariableDependency(id, needed_moose_vars, tid);
2506  _residual_materials.updateBoundaryMatPropDependency(id, needed_mat_props, tid);
2507  }
2508  }
2509 
2510  const std::set<MooseVariableFEBase *> & current_active_elemental_moose_variables =
2512  needed_moose_vars.insert(current_active_elemental_moose_variables.begin(),
2513  current_active_elemental_moose_variables.end());
2514 
2515  const std::set<unsigned int> & current_active_material_properties =
2517  needed_mat_props.insert(current_active_material_properties.begin(),
2518  current_active_material_properties.end());
2519 
2520  setActiveElementalMooseVariables(needed_moose_vars, tid);
2521  setActiveMaterialProperties(needed_mat_props, tid);
2522 }
2523 
2524 void
2525 FEProblemBase::reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful)
2526 {
2527  if (hasActiveMaterialProperties(tid))
2528  {
2529  auto && elem = _assembly[tid]->elem();
2530  unsigned int n_points = _assembly[tid]->qRule()->n_points();
2531  _material_data[tid]->resize(n_points);
2532 
2533  // Only swap if requested
2534  if (swap_stateful)
2535  _material_data[tid]->swap(*elem);
2536 
2537  if (_discrete_materials.hasActiveBlockObjects(blk_id, tid))
2539 
2541  _material_data[tid]->reinit(_jacobian_materials.getActiveBlockObjects(blk_id, tid));
2542 
2544  _material_data[tid]->reinit(_residual_materials.getActiveBlockObjects(blk_id, tid));
2545  }
2546 }
2547 
2548 void
2550 {
2551  if (hasActiveMaterialProperties(tid))
2552  {
2553  auto && elem = _assembly[tid]->elem();
2554  unsigned int side = _assembly[tid]->side();
2555  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2556 
2557  _bnd_material_data[tid]->resize(n_points);
2558 
2559  if (swap_stateful && !_bnd_material_data[tid]->isSwapped())
2560  _bnd_material_data[tid]->swap(*elem, side);
2561 
2562  if (_discrete_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2563  _bnd_material_data[tid]->reset(
2564  _discrete_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2565 
2566  if (_jacobian_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
2568  _bnd_material_data[tid]->reinit(
2569  _jacobian_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2570 
2571  if (_residual_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
2573  _bnd_material_data[tid]->reinit(
2574  _residual_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2575  }
2576 }
2577 
2578 void
2580 {
2581  if (hasActiveMaterialProperties(tid))
2582  {
2583  // NOTE: this will not work with h-adaptivity
2584  auto && neighbor = _assembly[tid]->neighbor();
2585  unsigned int neighbor_side = neighbor->which_neighbor_am_i(_assembly[tid]->elem());
2586  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2587  _neighbor_material_data[tid]->resize(n_points);
2588 
2589  // Only swap if requested
2590  if (swap_stateful)
2591  _neighbor_material_data[tid]->swap(*neighbor, neighbor_side);
2592 
2593  if (_discrete_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
2594  _neighbor_material_data[tid]->reset(
2595  _discrete_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2596 
2597  if (_jacobian_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
2599  _neighbor_material_data[tid]->reinit(
2600  _jacobian_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2601 
2602  if (_residual_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
2604  _neighbor_material_data[tid]->reinit(
2605  _residual_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));
2606  }
2607 }
2608 
2609 void
2610 FEProblemBase::reinitMaterialsBoundary(BoundaryID boundary_id, THREAD_ID tid, bool swap_stateful)
2611 {
2612  if (hasActiveMaterialProperties(tid))
2613  {
2614  auto && elem = _assembly[tid]->elem();
2615  unsigned int side = _assembly[tid]->side();
2616  unsigned int n_points = _assembly[tid]->qRuleFace()->n_points();
2617  _bnd_material_data[tid]->resize(n_points);
2618 
2619  if (swap_stateful && !_bnd_material_data[tid]->isSwapped())
2620  _bnd_material_data[tid]->swap(*elem, side);
2621 
2622  if (_discrete_materials.hasActiveBoundaryObjects(boundary_id, tid))
2623  _bnd_material_data[tid]->reset(
2624  _discrete_materials.getActiveBoundaryObjects(boundary_id, tid));
2625 
2626  if (_jacobian_materials.hasActiveBoundaryObjects(boundary_id, tid) &&
2628  _bnd_material_data[tid]->reinit(
2629  _jacobian_materials.getActiveBoundaryObjects(boundary_id, tid));
2630 
2631  if (_residual_materials.hasActiveBoundaryObjects(boundary_id, tid) &&
2633  _bnd_material_data[tid]->reinit(
2634  _residual_materials.getActiveBoundaryObjects(boundary_id, tid));
2635  }
2636 }
2637 
2638 void
2640 {
2641  auto && elem = _assembly[tid]->elem();
2642  _material_data[tid]->swapBack(*elem);
2643 }
2644 
2645 void
2647 {
2648  auto && elem = _assembly[tid]->elem();
2649  unsigned int side = _assembly[tid]->side();
2650  _bnd_material_data[tid]->swapBack(*elem, side);
2651 }
2652 
2653 void
2655 {
2656  // NOTE: this will not work with h-adaptivity
2657  auto && neighbor = _assembly[tid]->neighbor();
2658  unsigned int neighbor_side = neighbor->which_neighbor_am_i(_assembly[tid]->elem());
2659  _neighbor_material_data[tid]->swapBack(*neighbor, neighbor_side);
2660 }
2661 
2662 void
2663 FEProblemBase::initPostprocessorData(const std::string & name)
2664 {
2665  _pps_data.init(name);
2666 }
2667 
2668 void
2670 {
2671  _vpps_data.init(name);
2672 }
2673 
2674 void
2676  const std::string & name,
2677  InputParameters parameters)
2678 {
2679  // Check for name collision
2680  if (hasUserObject(name))
2681  mooseError(std::string("A UserObject with the name \"") + name +
2682  "\" already exists. You may not add a Postprocessor by the same name.");
2683 
2684  addUserObject(pp_name, name, parameters);
2686 }
2687 
2688 void
2690  const std::string & name,
2691  InputParameters parameters)
2692 {
2693  // Check for name collision
2694  if (hasUserObject(name))
2695  mooseError(std::string("A UserObject with the name \"") + name +
2696  "\" already exists. You may not add a VectorPostprocessor by the same name.");
2697 
2698  addUserObject(pp_name, name, parameters);
2700 }
2701 
2702 void
2703 FEProblemBase::addUserObject(std::string user_object_name,
2704  const std::string & name,
2705  InputParameters parameters)
2706 {
2707  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2708  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
2709  else
2710  {
2711  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
2712  {
2713  // We allow UserObjects to request that they use_displaced_mesh,
2714  // but then be overridden when no displacements variables are
2715  // provided in the Mesh block. If that happened, update the value
2716  // of use_displaced_mesh appropriately for this UserObject.
2717  if (parameters.have_parameter<bool>("use_displaced_mesh"))
2718  parameters.set<bool>("use_displaced_mesh") = false;
2719  }
2720 
2721  parameters.set<SubProblem *>("_subproblem") = this;
2722  }
2723 
2724  UserObject * primary = nullptr;
2725  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2726  {
2727  // Create the UserObject
2728  std::shared_ptr<UserObject> user_object =
2729  _factory.create<UserObject>(user_object_name, name, parameters, tid);
2730  if (tid == 0)
2731  primary = user_object.get();
2732  else
2733  user_object->setPrimaryThreadCopy(primary);
2734 
2735  // TODO: delete this line after apps have been updated to not call getUserObjects
2736  _all_user_objects.addObject(user_object, tid);
2737 
2738  theWarehouse().add(user_object, "UserObject");
2739 
2740  // Attempt to create all the possible UserObject types
2741  auto euo = std::dynamic_pointer_cast<ElementUserObject>(user_object);
2742  auto suo = std::dynamic_pointer_cast<SideUserObject>(user_object);
2743  auto isuo = std::dynamic_pointer_cast<InternalSideUserObject>(user_object);
2744  auto iuo = std::dynamic_pointer_cast<InterfaceUserObject>(user_object);
2745  auto nuo = std::dynamic_pointer_cast<NodalUserObject>(user_object);
2746  auto guo = std::dynamic_pointer_cast<GeneralUserObject>(user_object);
2747  auto tguo = std::dynamic_pointer_cast<ThreadedGeneralUserObject>(user_object);
2748 
2749  // Account for displaced mesh use
2750  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
2751  {
2752  if (euo || nuo)
2753  _reinit_displaced_elem = true;
2754  else if (suo || iuo)
2755  // shouldn't we add isuo
2756  _reinit_displaced_face = true;
2757  }
2758 
2759  if (guo && !tguo)
2760  break;
2761  }
2762 }
2763 
2764 const UserObject &
2765 FEProblemBase::getUserObjectBase(const std::string & name) const
2766 {
2767  std::vector<UserObject *> objs;
2768  theWarehouse().query().condition<AttribThread>(0).condition<AttribName>(name).queryInto(objs);
2769  if (objs.empty())
2770  mooseError("Unable to find user object with name '" + name + "'");
2771  return *(objs[0]);
2772 }
2773 
2774 bool
2775 FEProblemBase::hasUserObject(const std::string & name) const
2776 {
2777  std::vector<UserObject *> objs;
2778  theWarehouse().query().condition<AttribThread>(0).condition<AttribName>(name).queryInto(objs);
2779  return !objs.empty();
2780 }
2781 
2782 bool
2783 FEProblemBase::hasPostprocessor(const std::string & name)
2784 {
2786 }
2787 
2789 FEProblemBase::getPostprocessorValue(const PostprocessorName & name)
2790 {
2792 }
2793 
2796 {
2798 }
2799 
2802 {
2804 }
2805 
2808 {
2809  return _vpps_data;
2810 }
2811 
2812 bool
2813 FEProblemBase::hasVectorPostprocessor(const std::string & name)
2814 {
2816 }
2817 
2819 FEProblemBase::getVectorPostprocessorValue(const VectorPostprocessorName & name,
2820  const std::string & vector_name)
2821 {
2822  mooseDeprecated("getVectorPostprocessorValue() is DEPRECATED: Use the new version where you need "
2823  "to specify whether or not the vector must be broadcast");
2824 
2825  // The false means that we're not going to ask for this value to be broadcast
2826  // This mimics the old behavior - but is unsafe
2827  return _vpps_data.getVectorPostprocessorValue(name, vector_name, false);
2828 }
2829 
2832  const std::string & vector_name)
2833 {
2834  mooseDeprecated("getVectorPostprocessorValue() is DEPRECATED: Use the new version where you need "
2835  "to specify whether or not the vector must be broadcast");
2836 
2837  // The false means that we're not going to ask for this value to be broadcast
2838  // This mimics the old behavior - but is unsafe
2839  return _vpps_data.getVectorPostprocessorValueOld(name, vector_name, false);
2840 }
2841 
2843 FEProblemBase::getVectorPostprocessorValue(const VectorPostprocessorName & name,
2844  const std::string & vector_name,
2845  bool needs_broadcast)
2846 {
2847  return _vpps_data.getVectorPostprocessorValue(name, vector_name, needs_broadcast);
2848 }
2849 
2852  const std::string & vector_name,
2853  bool needs_broadcast)
2854 {
2855  return _vpps_data.getVectorPostprocessorValueOld(name, vector_name, needs_broadcast);
2856 }
2857 
2859 FEProblemBase::getScatterVectorPostprocessorValue(const VectorPostprocessorName & name,
2860  const std::string & vector_name)
2861 {
2862  return _vpps_data.getScatterVectorPostprocessorValue(name, vector_name);
2863 }
2864 
2866 FEProblemBase::getScatterVectorPostprocessorValueOld(const VectorPostprocessorName & name,
2867  const std::string & vector_name)
2868 {
2870 }
2871 
2873 FEProblemBase::declareVectorPostprocessorVector(const VectorPostprocessorName & name,
2874  const std::string & vector_name,
2875  bool contains_complete_history,
2876  bool is_broadcast)
2877 {
2878  return _vpps_data.declareVector(name, vector_name, contains_complete_history, is_broadcast);
2879 }
2880 
2881 const std::vector<std::pair<std::string, VectorPostprocessorData::VectorPostprocessorState>> &
2883 {
2884  return _vpps_data.vectors(vpp_name);
2885 }
2886 
2887 void
2889 {
2890  for (const auto & it : _multi_apps)
2891  {
2892  const auto & objects = it.second.getActiveObjects();
2893  for (const auto & obj : objects)
2894  obj->parentOutputPositionChanged();
2895  }
2896 }
2897 
2898 void
2900 {
2902  computeMarkers();
2903 }
2904 
2905 void
2907 {
2908  // Initialize indicator aux variable fields
2910  {
2911  TIME_SECTION(_compute_indicators_timer);
2912 
2913  std::vector<std::string> fields;
2914 
2915  // Indicator Fields
2916  const auto & indicators = _indicators.getActiveObjects();
2917  for (const auto & indicator : indicators)
2918  fields.push_back(indicator->name());
2919 
2920  // InternalSideIndicator Fields
2921  const auto & internal_indicators = _internal_side_indicators.getActiveObjects();
2922  for (const auto & internal_indicator : internal_indicators)
2923  fields.push_back(internal_indicator->name());
2924 
2925  _aux->zeroVariables(fields);
2926 
2927  // compute Indicators
2928  ComputeIndicatorThread cit(*this);
2929  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cit);
2930  _aux->solution().close();
2931  _aux->update();
2932 
2933  ComputeIndicatorThread finalize_cit(*this, true);
2934  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), finalize_cit);
2935  _aux->solution().close();
2936  _aux->update();
2937  }
2938 }
2939 
2940 void
2942 {
2943  if (_markers.hasActiveObjects())
2944  {
2945  TIME_SECTION(_compute_markers_timer);
2946 
2947  std::vector<std::string> fields;
2948 
2949  // Marker Fields
2950  const auto & markers = _markers.getActiveObjects();
2951  for (const auto & marker : markers)
2952  fields.push_back(marker->name());
2953 
2954  _aux->zeroVariables(fields);
2955 
2957 
2958  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
2959  {
2960  const auto & markers = _markers.getActiveObjects(tid);
2961  for (const auto & marker : markers)
2962  marker->markerSetup();
2963  }
2964 
2965  ComputeMarkerThread cmt(*this);
2966  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cmt);
2967 
2968  _aux->solution().close();
2969  _aux->update();
2970  }
2971 }
2972 
2973 const ExecFlagType &
2975 {
2976  return _current_execute_on_flag;
2977 }
2978 
2979 void
2981 {
2982  _current_execute_on_flag = flag;
2983 }
2984 
2985 void
2987 {
2988  // Set the current flag
2989  setCurrentExecuteOnFlag(exec_type);
2990  if (exec_type == EXEC_NONLINEAR)
2991  {
2993  if (_displaced_problem)
2994  _displaced_problem->setCurrentlyComputingJacobian(true);
2995  }
2996 
2997  // Samplers
2998  if (exec_type != EXEC_INITIAL)
2999  executeSamplers(exec_type);
3000 
3001  // Pre-aux UserObjects
3002  computeUserObjects(exec_type, Moose::PRE_AUX);
3003 
3004  // AuxKernels
3005  computeAuxiliaryKernels(exec_type);
3006 
3007  // Post-aux UserObjects
3008  computeUserObjects(exec_type, Moose::POST_AUX);
3009 
3010  // Controls
3011  if (exec_type != EXEC_INITIAL)
3012  executeControls(exec_type);
3013 
3014  // Return the current flag to None
3017  if (_displaced_problem)
3018  _displaced_problem->setCurrentlyComputingJacobian(false);
3019 }
3020 
3021 void
3023 {
3024  _aux->compute(type);
3025 }
3026 
3027 // Finalize, threadJoin, and update PP values of Elemental/Nodal/Side/InternalSideUserObjects
3028 void
3030 {
3031  std::vector<UserObject *> objs;
3032  query.queryInto(objs);
3033  if (!isgen)
3034  {
3035  // join all threaded user objects (i.e. not regular general user objects) to the primary thread
3036  for (auto obj : objs)
3037  if (obj->primaryThreadCopy())
3038  obj->primaryThreadCopy()->threadJoin(*obj);
3039  }
3040 
3041  query.condition<AttribThread>(0).queryInto(objs);
3042 
3043  // finalize objects and retrieve/store any postproessor values
3044  for (auto obj : objs)
3045  {
3046  if (isgen && dynamic_cast<ThreadedGeneralUserObject *>(obj))
3047  continue;
3048  if (isgen)
3049  {
3050  // general user objects are not run in their own threaded loop object - so run them here
3051  obj->initialize();
3052  obj->execute();
3053  }
3054 
3055  obj->finalize();
3056 
3057  // These have to be stored piecemeal (with every call to this function) because general
3058  // postprocessors (which run last after other userobjects have been completed) might depend on
3059  // them being stored. This wouldn't be a problem if all userobjects satisfied the dependency
3060  // resolver interface and could be sorted appropriately with the general userobjects, but they
3061  // don't.
3062  auto pp = dynamic_cast<Postprocessor *>(obj);
3063  if (pp)
3064  _pps_data.storeValue(pp->PPName(), pp->getValue());
3065  auto vpp = dynamic_cast<VectorPostprocessor *>(obj);
3066  if (vpp)
3067  _vpps_data.broadcastScatterVectors(vpp->PPName());
3068  }
3069 }
3070 
3071 void
3073 {
3075  .query()
3076  .condition<AttribSystem>("UserObject")
3077  .condition<AttribExecOns>(type)
3078  .condition<AttribName>(name);
3080 }
3081 
3082 void
3084 {
3085  TheWarehouse::Query query =
3086  theWarehouse().query().condition<AttribSystem>("UserObject").condition<AttribExecOns>(type);
3087  computeUserObjectsInternal(type, group, query);
3088 }
3089 
3090 void
3092  const Moose::AuxGroup & group,
3093  TheWarehouse::Query & query)
3094 {
3095  if (group == Moose::PRE_IC)
3096  query.condition<AttribPreIC>(true);
3097  else if (group == Moose::PRE_AUX)
3098  query.condition<AttribPreAux>(true);
3099  else if (group == Moose::POST_AUX)
3100  query.condition<AttribPreAux>(false);
3101 
3102  std::vector<GeneralUserObject *> genobjs;
3103  query.clone().condition<AttribInterfaces>(Interfaces::GeneralUserObject).queryInto(genobjs);
3104 
3105  std::vector<UserObject *> userobjs;
3106  query.clone()
3110  .queryInto(userobjs);
3111 
3112  std::vector<UserObject *> tgobjs;
3113  query.clone()
3115  .queryInto(tgobjs);
3116 
3117  std::vector<UserObject *> nodal;
3118  query.clone().condition<AttribInterfaces>(Interfaces::NodalUserObject).queryInto(nodal);
3119 
3120  if (userobjs.empty() && genobjs.empty() && tgobjs.empty() && nodal.empty())
3121  return;
3122 
3123  TIME_SECTION(_compute_user_objects_timer);
3124 
3125  // Start the timer here since we have at least one active user object
3126  std::string compute_uo_tag = "computeUserObjects(" + Moose::stringify(type) + ")";
3127 
3128  // Perform Residual/Jacobian setups
3129  if (type == EXEC_LINEAR)
3130  {
3131  for (auto obj : userobjs)
3132  obj->residualSetup();
3133  for (auto obj : nodal)
3134  obj->residualSetup();
3135  for (auto obj : tgobjs)
3136  obj->residualSetup();
3137  for (auto obj : genobjs)
3138  obj->residualSetup();
3139  }
3140  else if (type == EXEC_NONLINEAR)
3141  {
3142  for (auto obj : userobjs)
3143  obj->jacobianSetup();
3144  for (auto obj : nodal)
3145  obj->jacobianSetup();
3146  for (auto obj : tgobjs)
3147  obj->jacobianSetup();
3148  for (auto obj : genobjs)
3149  obj->jacobianSetup();
3150  }
3151 
3152  for (auto obj : userobjs)
3153  obj->initialize();
3154 
3155  // Execute Elemental/Side/InternalSideUserObjects
3156  if (!userobjs.empty())
3157  {
3158  // non-nodal user objects have to be run separately before the nodal user objects run
3159  // because some nodal user objects (NodalNormal related) depend on elemental user objects :-(
3160  ComputeUserObjectsThread cppt(*this, getNonlinearSystemBase(), query);
3161  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), cppt);
3162 
3163  // There is one instance in rattlesnake where an elemental user object's finalize depends
3164  // on a side user object having been finalized first :-(
3169  }
3170 
3171  // Execute NodalUserObjects
3172  // BISON has an axial reloc elemental user object that has a finalize func that depends on a
3173  // nodal user object's prev value. So we can't initialize this until after elemental objects
3174  // have been finalized :-(
3175  for (auto obj : nodal)
3176  obj->initialize();
3177  if (query.clone().condition<AttribInterfaces>(Interfaces::NodalUserObject).count() > 0)
3178  {
3179  ComputeNodalUserObjectsThread cnppt(*this, query);
3180  Threads::parallel_reduce(*_mesh.getLocalNodeRange(), cnppt);
3182  }
3183 
3184  for (auto obj : tgobjs)
3185  obj->initialize();
3186  std::vector<GeneralUserObject *> tguos_zero;
3187  query.clone()
3188  .condition<AttribThread>(0)
3189  .condition<AttribInterfaces>(Interfaces::ThreadedGeneralUserObject)
3190  .queryInto(tguos_zero);
3191  for (auto obj : tguos_zero)
3192  {
3193  std::vector<GeneralUserObject *> tguos;
3194  auto q = query.clone()
3195  .condition<AttribName>(obj->name())
3196  .condition<AttribInterfaces>(Interfaces::ThreadedGeneralUserObject);
3197  q.queryInto(tguos);
3198 
3200  Threads::parallel_reduce(GeneralUserObjectRange(tguos.begin(), tguos.end()), ctguot);
3201  joinAndFinalize(q);
3202  }
3203 
3205 }
3206 
3207 void
3209 {
3210  if (_control_warehouse[exec_type].hasActiveObjects())
3211  {
3212  TIME_SECTION(_execute_controls_timer);
3213 
3215 
3216  auto controls_wh = _control_warehouse[exec_type];
3217  // Add all of the dependencies into the resolver and sort them
3218  for (const auto & it : controls_wh.getActiveObjects())
3219  {
3220  // Make sure an item with no dependencies comes out too!
3221  resolver.addItem(it);
3222 
3223  std::vector<std::string> & dependent_controls = it->getDependencies();
3224  for (const auto & depend_name : dependent_controls)
3225  {
3226  if (controls_wh.hasActiveObject(depend_name))
3227  {
3228  auto dep_control = controls_wh.getActiveObject(depend_name);
3229  resolver.insertDependency(it, dep_control);
3230  }
3231  else
3232  mooseError("The Control \"",
3233  depend_name,
3234  "\" was not created, did you make a "
3235  "spelling mistake or forget to include it "
3236  "in your input file?");
3237  }
3238  }
3239 
3240  const auto & ordered_controls = resolver.getSortedValues();
3241 
3242  if (!ordered_controls.empty())
3243  {
3244  _control_warehouse.setup(exec_type);
3245  // Run the controls in the proper order
3246  for (const auto & control : ordered_controls)
3247  control->execute();
3248  }
3249  }
3250 }
3251 
3252 void
3254 {
3255  // TODO: This should be done in a threaded loop, but this should be super quick so for now
3256  // do a serial loop.
3257  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3258  {
3259  const auto & objects = _samplers[exec_type].getActiveObjects(tid);
3260  if (!objects.empty())
3261  {
3262  TIME_SECTION(_execute_samplers_timer);
3263 
3264  _samplers.setup(exec_type);
3265  for (auto & sampler : objects)
3266  sampler->execute();
3267  }
3268  }
3269 }
3270 
3271 void
3273 {
3274  TIME_SECTION(_update_active_objects_timer);
3275 
3276  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
3277  {
3278  _nl->updateActive(tid);
3279  _aux->updateActive(tid);
3282  _markers.updateActive(tid);
3287  _samplers.updateActive(tid);
3288  }
3289 
3296 }
3297 
3298 void
3300 {
3301  //<< "Object " << a->name() << " -> " << b->name() << std::endl;
3302 }
3303 
3304 void
3306 {
3308 
3309  // Need to see if _any_ processor has ghosted elems or geometry objects.
3310  bool needs_reinit = !_ghosted_elems.empty();
3311  needs_reinit = needs_reinit || !_geometric_search_data._nearest_node_locators.empty() ||
3313  needs_reinit =
3314  needs_reinit ||
3315  (_displaced_problem && !_displaced_problem->geomSearchData()._nearest_node_locators.empty());
3316  _communicator.max(needs_reinit);
3317 
3318  if (needs_reinit)
3319  {
3320  // Call reinit to get the ghosted vectors correct now that some geometric search has been done
3321  _eq.reinit();
3322 
3323  if (_displaced_mesh)
3324  _displaced_problem->es().reinit();
3325  }
3326 }
3327 
3328 void
3329 FEProblemBase::addDamper(std::string damper_name,
3330  const std::string & name,
3331  InputParameters parameters)
3332 {
3333  parameters.set<SubProblem *>("_subproblem") = this;
3334  parameters.set<SystemBase *>("_sys") = _nl.get();
3335 
3336  _has_dampers = true;
3337  _nl->addDamper(damper_name, name, parameters);
3338 }
3339 
3340 void
3342 {
3343  _nl->setupDampers();
3344 }
3345 
3346 void
3347 FEProblemBase::addIndicator(std::string indicator_name,
3348  const std::string & name,
3349  InputParameters parameters)
3350 {
3351  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
3352  {
3353  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3354  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3355  _reinit_displaced_elem = true;
3356  }
3357  else
3358  {
3359  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
3360  {
3361  // We allow Indicators to request that they use_displaced_mesh,
3362  // but then be overridden when no displacements variables are
3363  // provided in the Mesh block. If that happened, update the value
3364  // of use_displaced_mesh appropriately for this Indicator.
3365  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3366  parameters.set<bool>("use_displaced_mesh") = false;
3367  }
3368 
3369  parameters.set<SubProblem *>("_subproblem") = this;
3370  parameters.set<SystemBase *>("_sys") = _aux.get();
3371  }
3372 
3373  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
3374  {
3375  std::shared_ptr<Indicator> indicator =
3376  _factory.create<Indicator>(indicator_name, name, parameters, tid);
3377 
3378  std::shared_ptr<InternalSideIndicator> isi =
3380  if (isi)
3382  else
3383  _indicators.addObject(indicator, tid);
3384  }
3385 }
3386 
3387 void
3388 FEProblemBase::addMarker(std::string marker_name,
3389  const std::string & name,
3390  InputParameters parameters)
3391 {
3392  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
3393  {
3394  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3395  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3396  _reinit_displaced_elem = true;
3397  }
3398  else
3399  {
3400  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
3401  {
3402  // We allow Markers to request that they use_displaced_mesh,
3403  // but then be overridden when no displacements variables are
3404  // provided in the Mesh block. If that happened, update the value
3405  // of use_displaced_mesh appropriately for this Marker.
3406  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3407  parameters.set<bool>("use_displaced_mesh") = false;
3408  }
3409 
3410  parameters.set<SubProblem *>("_subproblem") = this;
3411  parameters.set<SystemBase *>("_sys") = _aux.get();
3412  }
3413 
3414  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
3415  {
3416  std::shared_ptr<Marker> marker = _factory.create<Marker>(marker_name, name, parameters, tid);
3417  _markers.addObject(marker, tid);
3418  }
3419 }
3420 
3421 void
3422 FEProblemBase::addMultiApp(const std::string & multi_app_name,
3423  const std::string & name,
3424  InputParameters parameters)
3425 {
3426  parameters.set<MPI_Comm>("_mpi_comm") = _communicator.get();
3427  parameters.set<std::shared_ptr<CommandLine>>("_command_line") = _app.commandLine();
3428 
3429  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
3430  {
3431  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3432  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3433  _reinit_displaced_elem = true;
3434  }
3435  else
3436  {
3437  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
3438  {
3439  // We allow MultiApps to request that they use_displaced_mesh,
3440  // but then be overridden when no displacements variables are
3441  // provided in the Mesh block. If that happened, update the value
3442  // of use_displaced_mesh appropriately for this MultiApp.
3443  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3444  parameters.set<bool>("use_displaced_mesh") = false;
3445  }
3446 
3447  parameters.set<SubProblem *>("_subproblem") = this;
3448  parameters.set<SystemBase *>("_sys") = _aux.get();
3449  }
3450 
3451  std::shared_ptr<MultiApp> multi_app = _factory.create<MultiApp>(multi_app_name, name, parameters);
3452 
3453  multi_app->setupPositions();
3454 
3455  _multi_apps.addObject(multi_app);
3456 
3457  // Store TranseintMultiApp objects in another container, this is needed for calling computeDT
3458  std::shared_ptr<TransientMultiApp> trans_multi_app =
3460  if (trans_multi_app)
3461  _transient_multi_apps.addObject(trans_multi_app);
3462 }
3463 
3464 bool
3466 {
3467  return _multi_apps[type].hasActiveObjects();
3468 }
3469 
3470 bool
3471 FEProblemBase::hasMultiApp(const std::string & multi_app_name) const
3472 {
3473  return _multi_apps.hasActiveObject(multi_app_name);
3474 }
3475 
3476 std::shared_ptr<MultiApp>
3477 FEProblemBase::getMultiApp(const std::string & multi_app_name) const
3478 {
3479  return _multi_apps.getObject(multi_app_name);
3480 }
3481 
3482 void
3484 {
3485  bool to_multiapp = direction == MultiAppTransfer::TO_MULTIAPP;
3486  std::string string_direction = to_multiapp ? " To " : " From ";
3487  const MooseObjectWarehouse<Transfer> & wh =
3489 
3490  if (wh.hasActiveObjects())
3491  {
3492  TIME_SECTION(_exec_multi_app_transfers_timer);
3493 
3494  const auto & transfers = wh.getActiveObjects();
3495 
3496  _console << COLOR_CYAN << "\nStarting Transfers on " << Moose::stringify(type)
3497  << string_direction << "MultiApps" << COLOR_DEFAULT << std::endl;
3498  for (const auto & transfer : transfers)
3499  transfer->execute();
3500 
3502 
3503  _console << COLOR_CYAN << "Transfers on " << Moose::stringify(type) << " Are Finished\n"
3504  << COLOR_DEFAULT << std::endl;
3505  }
3506  else if (_multi_apps[type].getActiveObjects().size())
3507  _console << COLOR_CYAN << "\nNo Transfers on " << Moose::stringify(type) << " To MultiApps\n"
3508  << COLOR_DEFAULT << std::endl;
3509 }
3510 
3511 std::vector<std::shared_ptr<Transfer>>
3513 {
3517  return wh.getActiveObjects();
3518 }
3519 
3522 {
3523  if (direction == MultiAppTransfer::TO_MULTIAPP)
3524  return _to_multi_app_transfers;
3525  else
3527 }
3528 
3529 bool
3531 {
3532  // Active MultiApps
3533  const std::vector<MooseSharedPointer<MultiApp>> & multi_apps =
3535 
3536  // Do anything that needs to be done to Apps before transfers
3537  for (const auto & multi_app : multi_apps)
3538  multi_app->preTransfer(_dt, _time);
3539 
3540  // Execute Transfers _to_ MultiApps
3542 
3543  // Execute MultiApps
3544  if (multi_apps.size())
3545  {
3546  TIME_SECTION(_exec_multi_apps_timer);
3547 
3548  _console << COLOR_CYAN << "\nExecuting MultiApps on " << Moose::stringify(type) << COLOR_DEFAULT
3549  << std::endl;
3550 
3551  bool success = true;
3552 
3553  for (const auto & multi_app : multi_apps)
3554  {
3555  success = multi_app->solveStep(_dt, _time, auto_advance);
3556  // no need to finish executing the subapps if one fails
3557  if (!success)
3558  break;
3559  }
3560 
3562 
3563  _communicator.min(success);
3564 
3565  if (!success)
3566  return false;
3567 
3568  _console << COLOR_CYAN << "Finished Executing MultiApps on " << Moose::stringify(type) << "\n"
3569  << COLOR_DEFAULT << std::endl;
3570  }
3571 
3572  // Execute Transfers _from_ MultiApps
3574 
3575  // If we made it here then everything passed
3576  return true;
3577 }
3578 
3579 void
3581 {
3582  const auto & multi_apps = _multi_apps.getActiveObjects();
3583 
3584  for (const auto & multi_app : multi_apps)
3585  multi_app->finalize();
3586 }
3587 
3588 void
3590 {
3591  const auto & multi_apps = _multi_apps.getActiveObjects();
3592 
3593  for (const auto & multi_app : multi_apps)
3594  multi_app->postExecute();
3595 }
3596 
3597 void
3599 {
3600  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3601 
3602  if (multi_apps.size())
3603  for (const auto & multi_app : multi_apps)
3604  multi_app->incrementTStep(_time);
3605 }
3606 
3607 void
3609 {
3610  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3611 
3612  if (multi_apps.size())
3613  {
3614  _console << COLOR_CYAN << "\nAdvancing MultiApps on " << type.name() << COLOR_DEFAULT
3615  << std::endl;
3616 
3617  for (const auto & multi_app : multi_apps)
3618  multi_app->finishStep();
3619 
3621 
3622  _console << COLOR_CYAN << "Finished Advancing MultiApps on " << type.name() << "\n"
3623  << COLOR_DEFAULT << std::endl;
3624  }
3625 }
3626 
3627 void
3629 {
3630  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3631 
3632  if (multi_apps.size())
3633  {
3634  TIME_SECTION(_backup_multi_apps_timer);
3635 
3636  _console << COLOR_CYAN << "\nBacking Up MultiApps on " << type.name() << COLOR_DEFAULT
3637  << std::endl;
3638 
3639  for (const auto & multi_app : multi_apps)
3640  multi_app->backup();
3641 
3643 
3644  _console << COLOR_CYAN << "Finished Backing Up MultiApps on " << type.name() << "\n"
3645  << COLOR_DEFAULT << std::endl;
3646  }
3647 }
3648 
3649 void
3651 {
3652  const auto & multi_apps = _multi_apps[type].getActiveObjects();
3653 
3654  if (multi_apps.size())
3655  {
3656  if (force)
3657  _console << COLOR_CYAN << "\nRestoring Multiapps on " << type.name()
3658  << " because of solve failure!" << COLOR_DEFAULT << std::endl;
3659  else
3660  _console << COLOR_CYAN << "\nRestoring MultiApps on " << type.name() << COLOR_DEFAULT
3661  << std::endl;
3662 
3663  for (const auto & multi_app : multi_apps)
3664  if (force || multi_app->needsRestoration())
3665  multi_app->restore();
3666 
3668 
3669  _console << COLOR_CYAN << "Finished Restoring MultiApps on " << type.name() << "\n"
3670  << COLOR_DEFAULT << std::endl;
3671  }
3672 }
3673 
3674 Real
3676 {
3677  const auto & multi_apps = _transient_multi_apps[type].getActiveObjects();
3678 
3679  Real smallest_dt = std::numeric_limits<Real>::max();
3680 
3681  for (const auto & multi_app : multi_apps)
3682  smallest_dt = std::min(smallest_dt, multi_app->computeDT());
3683 
3684  return smallest_dt;
3685 }
3686 
3687 void
3689 {
3690  if (_transfers[type].hasActiveObjects())
3691  {
3692  const auto & transfers = _transfers[type].getActiveObjects();
3693  for (const auto & transfer : transfers)
3694  transfer->execute();
3695  }
3696 }
3697 
3698 void
3699 FEProblemBase::addTransfer(const std::string & transfer_name,
3700  const std::string & name,
3701  InputParameters parameters)
3702 {
3703  if (_displaced_problem && parameters.get<bool>("use_displaced_mesh"))
3704  {
3705  parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
3706  parameters.set<SystemBase *>("_sys") = &_displaced_problem->auxSys();
3707  _reinit_displaced_elem = true;
3708  }
3709  else
3710  {
3711  if (_displaced_problem == nullptr && parameters.get<bool>("use_displaced_mesh"))
3712  {
3713  // We allow Transfers to request that they use_displaced_mesh,
3714  // but then be overridden when no displacements variables are
3715  // provided in the Mesh block. If that happened, update the value
3716  // of use_displaced_mesh appropriately for this Transfer.
3717  if (parameters.have_parameter<bool>("use_displaced_mesh"))
3718  parameters.set<bool>("use_displaced_mesh") = false;
3719  }
3720 
3721  parameters.set<SubProblem *>("_subproblem") = this;
3722  parameters.set<SystemBase *>("_sys") = _aux.get();
3723  }
3724 
3725  // Handle the "SAME_AS_MULTIAPP" execute option. The get method is used to test for the
3726  // flag so the set by user flag is not reset, calling set with the true flag causes the set
3727  // by user status to be reset, which should only be done if the EXEC_SAME_AS_MULTIAPP is
3728  // being applied to the object.
3729  if (parameters.get<ExecFlagEnum>("execute_on").contains(EXEC_SAME_AS_MULTIAPP))
3730  {
3731  ExecFlagEnum & exec_enum = parameters.set<ExecFlagEnum>("execute_on", true);
3732  std::shared_ptr<MultiApp> multiapp = getMultiApp(parameters.get<MultiAppName>("multi_app"));
3733  exec_enum = multiapp->getParam<ExecFlagEnum>("execute_on");
3734  }
3735 
3736  // Create the Transfer objects
3737  std::shared_ptr<Transfer> transfer = _factory.create<Transfer>(transfer_name, name, parameters);
3738 
3739  // Add MultiAppTransfer object
3740  std::shared_ptr<MultiAppTransfer> multi_app_transfer =
3742  if (multi_app_transfer)
3743  {
3744  if (multi_app_transfer->direction() == MultiAppTransfer::TO_MULTIAPP)
3745  _to_multi_app_transfers.addObject(multi_app_transfer);
3746  else
3747  _from_multi_app_transfers.addObject(multi_app_transfer);
3748  }
3749  else
3750  _transfers.addObject(transfer);
3751 }
3752 
3753 bool
3754 FEProblemBase::hasVariable(const std::string & var_name) const
3755 {
3756  if (_nl->hasVariable(var_name))
3757  return true;
3758  else if (_aux->hasVariable(var_name))
3759  return true;
3760  else
3761  return false;
3762 }
3763 
3766  const std::string & var_name,
3767  Moose::VarKindType expected_var_type,
3768  Moose::VarFieldType expected_var_field_type)
3769 {
3770  return getVariableHelper(tid, var_name, expected_var_type, expected_var_field_type, *_nl, *_aux);
3771 }
3772 
3773 MooseVariable &
3774 FEProblemBase::getStandardVariable(THREAD_ID tid, const std::string & var_name)
3775 {
3776  if (_nl->hasVariable(var_name))
3777  return _nl->getFieldVariable<Real>(tid, var_name);
3778  else if (!_aux->hasVariable(var_name))
3779  mooseError("Unknown variable " + var_name);
3780 
3781  return _aux->getFieldVariable<Real>(tid, var_name);
3782 }
3783 
3785 FEProblemBase::getVectorVariable(THREAD_ID tid, const std::string & var_name)
3786 {
3787  if (_nl->hasVariable(var_name))
3788  return _nl->getFieldVariable<RealVectorValue>(tid, var_name);
3789  else if (!_aux->hasVariable(var_name))
3790  mooseError("Unknown variable " + var_name);
3791 
3792  return _aux->getFieldVariable<RealVectorValue>(tid, var_name);
3793 }
3794 
3795 bool
3796 FEProblemBase::hasScalarVariable(const std::string & var_name) const
3797 {
3798  if (_nl->hasScalarVariable(var_name))
3799  return true;
3800  else if (_aux->hasScalarVariable(var_name))
3801  return true;
3802  else
3803  return false;
3804 }
3805 
3807 FEProblemBase::getScalarVariable(THREAD_ID tid, const std::string & var_name)
3808 {
3809  if (_nl->hasScalarVariable(var_name))
3810  return _nl->getScalarVariable(tid, var_name);
3811  else if (_aux->hasScalarVariable(var_name))
3812  return _aux->getScalarVariable(tid, var_name);
3813  else
3814  mooseError("Unknown variable " + var_name);
3815 }
3816 
3817 System &
3818 FEProblemBase::getSystem(const std::string & var_name)
3819 {
3820  if (_nl->hasVariable(var_name))
3821  return _nl->system();
3822  else if (_aux->hasVariable(var_name))
3823  return _aux->system();
3824  else
3825  mooseError("Unable to find a system containing the variable " + var_name);
3826 }
3827 
3828 void
3830 {
3832 
3833  if (_displaced_problem)
3834  _displaced_problem->setActiveFEVariableCoupleableMatrixTags(mtags, tid);
3835 }
3836 
3837 void
3839 {
3841 
3842  if (_displaced_problem)
3843  _displaced_problem->setActiveFEVariableCoupleableVectorTags(vtags, tid);
3844 }
3845 
3846 void
3848 {
3850 
3851  if (_displaced_problem)
3852  _displaced_problem->setActiveScalarVariableCoupleableMatrixTags(mtags, tid);
3853 }
3854 
3855 void
3857 {
3859 
3860  if (_displaced_problem)
3861  _displaced_problem->setActiveScalarVariableCoupleableVectorTags(vtags, tid);
3862 }
3863 
3864 void
3865 FEProblemBase::setActiveElementalMooseVariables(const std::set<MooseVariableFEBase *> & moose_vars,
3866  THREAD_ID tid)
3867 {
3869 
3870  if (_displaced_problem)
3871  _displaced_problem->setActiveElementalMooseVariables(moose_vars, tid);
3872 }
3873 
3874 void
3876 {
3878 
3879  if (_displaced_problem)
3880  _displaced_problem->clearActiveElementalMooseVariables(tid);
3881 }
3882 
3883 void
3885 {
3887 
3888  if (_displaced_problem)
3889  _displaced_problem->clearActiveFEVariableCoupleableMatrixTags(tid);
3890 }
3891 
3892 void
3894 {
3896 
3897  if (_displaced_problem)
3898  _displaced_problem->clearActiveFEVariableCoupleableVectorTags(tid);
3899 }
3900 
3901 void
3903 {
3905 
3906  if (_displaced_problem)
3907  _displaced_problem->clearActiveScalarVariableCoupleableMatrixTags(tid);
3908 }
3909 
3910 void
3912 {
3914 
3915  if (_displaced_problem)
3916  _displaced_problem->clearActiveScalarVariableCoupleableVectorTags(tid);
3917 }
3918 
3919 void
3920 FEProblemBase::setActiveMaterialProperties(const std::set<unsigned int> & mat_prop_ids,
3921  THREAD_ID tid)
3922 {
3923  SubProblem::setActiveMaterialProperties(mat_prop_ids, tid);
3924 
3925  if (_displaced_problem)
3926  _displaced_problem->setActiveMaterialProperties(mat_prop_ids, tid);
3927 }
3928 
3929 void
3931 {
3933 
3934  if (_displaced_problem)
3935  _displaced_problem->clearActiveMaterialProperties(tid);
3936 }
3937 
3938 void
3939 FEProblemBase::createQRules(QuadratureType type, Order order, Order volume_order, Order face_order)
3940 {
3941  if (order == INVALID_ORDER)
3942  {
3943  // automatically determine the integration order
3944  order = _nl->getMinQuadratureOrder();
3945  if (order < _aux->getMinQuadratureOrder())
3946  order = _aux->getMinQuadratureOrder();
3947  }
3948 
3949  if (volume_order == INVALID_ORDER)
3950  volume_order = order;
3951 
3952  if (face_order == INVALID_ORDER)
3953  face_order = order;
3954 
3955  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
3956  _assembly[tid]->createQRules(type, order, volume_order, face_order);
3957 
3958  if (_displaced_problem)
3959  _displaced_problem->createQRules(type, order, volume_order, face_order);
3960 
3961  // Find the maximum number of quadrature points
3962  {
3963  MaxQpsThread mqt(*this, type, std::max(order, volume_order), face_order);
3964  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mqt);
3965  _max_qps = mqt.max();
3967 
3968  // If we have more shape functions or more quadrature points on
3969  // another processor, then we may need to handle those elements
3970  // ourselves later after repartitioning.
3971  _communicator.max(_max_qps);
3972  _communicator.max(_max_shape_funcs);
3973  }
3974 
3975  unsigned int max_qpts = getMaxQps();
3976  for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
3977  {
3978  // the highest available order in libMesh is 43
3979  _scalar_zero[tid].resize(FORTYTHIRD, 0);
3980  _zero[tid].resize(max_qpts, 0);
3981  _ad_zero[tid].resize(max_qpts, 0);
3982  _grad_zero[tid].resize(max_qpts, RealGradient(0.));
3983  _ad_grad_zero[tid].resize(max_qpts, DualRealGradient(0));
3984  _second_zero[tid].resize(max_qpts, RealTensor(0.));
3985  _ad_second_zero[tid].resize(max_qpts, DualRealTensorValue(0));
3986  _second_phi_zero[tid].resize(max_qpts,
3987  std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));
3988  _vector_zero[tid].resize(max_qpts, RealGradient(0.));
3989  _vector_curl_zero[tid].resize(max_qpts, RealGradient(0.));
3990  }
3991 }
3992 
3993 void
3995 {
3996  _coupling = type;
3997 }
3998 
3999 void
4001 {
4002  // TODO: Deprecate method
4003 
4005  _cm.reset(cm);
4006 }
4007 
4008 void
4009 FEProblemBase::setCouplingMatrix(std::unique_ptr<CouplingMatrix> cm)
4010 {
4012  _cm = std::move(cm);
4013 }
4014 
4015 void
4017 {
4018  unsigned int n_vars = _nl->nVariables();
4019  _nonlocal_cm.resize(n_vars);
4020  const auto & vars = _nl->getVariables(0);
4021  const auto & nonlocal_kernel = _nonlocal_kernels.getObjects();
4022  const auto & nonlocal_integrated_bc = _nonlocal_integrated_bcs.getObjects();
4023  for (const auto & ivar : vars)
4024  {
4025  for (const auto & kernel : nonlocal_kernel)
4026  {
4027  unsigned int i = ivar->number();
4028  if (i == kernel->variable().number())
4029  for (const auto & jvar : vars)
4030  {
4031  const auto it = _var_dof_map.find(jvar->name());
4032  if (it != _var_dof_map.end())
4033  {
4034  unsigned int j = jvar->number();
4035  _nonlocal_cm(i, j) = 1;
4036  }
4037  }
4038  }
4039  for (const auto & integrated_bc : nonlocal_integrated_bc)
4040  {
4041  unsigned int i = ivar->number();
4042  if (i == integrated_bc->variable().number())
4043  for (const auto & jvar : vars)
4044  {
4045  const auto it = _var_dof_map.find(jvar->name());
4046  if (it != _var_dof_map.end())
4047  {
4048  unsigned int j = jvar->number();
4049  _nonlocal_cm(i, j) = 1;
4050  }
4051  }
4052  }
4053  }
4054 }
4055 
4056 bool
4057 FEProblemBase::areCoupled(unsigned int ivar, unsigned int jvar)
4058 {
4059  return (*_cm)(ivar, jvar);
4060 }
4061 
4062 std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> &
4064 {
4065  return _assembly[tid]->couplingEntries();
4066 }
4067 
4068 std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> &
4070 {
4071  return _assembly[tid]->nonlocalCouplingEntries();
4072 }
4073 
4074 void
4076 {
4077  if (_initialized)
4078  return;
4079 
4080  TIME_SECTION(_init_timer);
4081 
4082  unsigned int n_vars = _nl->nVariables();
4083  switch (_coupling)
4084  {
4085  case Moose::COUPLING_DIAG:
4086  _cm = libmesh_make_unique<CouplingMatrix>(n_vars);
4087  for (unsigned int i = 0; i < n_vars; i++)
4088  for (unsigned int j = 0; j < n_vars; j++)
4089  (*_cm)(i, j) = (i == j ? 1 : 0);
4090  break;
4091 
4092  // for full jacobian
4093  case Moose::COUPLING_FULL:
4094  _cm = libmesh_make_unique<CouplingMatrix>(n_vars);
4095  for (unsigned int i = 0; i < n_vars; i++)
4096  for (unsigned int j = 0; j < n_vars; j++)
4097  (*_cm)(i, j) = 1;
4098  break;
4099 
4101  // do nothing, _cm was already set through couplingMatrix() call
4102  break;
4103  }
4104 
4105  _nl->dofMap()._dof_coupling = _cm.get();
4106 
4107  // If there are no variables, make sure to pass a nullptr coupling
4108  // matrix, to avoid warnings about non-nullptr yet empty
4109  // CouplingMatrices.
4110  if (n_vars == 0)
4111  _nl->dofMap()._dof_coupling = nullptr;
4112 
4113  _nl->dofMap().attach_extra_sparsity_function(&extraSparsity, _nl.get());
4114  _nl->dofMap().attach_extra_send_list_function(&extraSendList, _nl.get());
4115  _aux->dofMap().attach_extra_send_list_function(&extraSendList, _aux.get());
4116 
4117  if (!_skip_nl_system_check && _solve && n_vars == 0)
4118  mooseError("No variables specified in the FEProblemBase '", name(), "'.");
4119 
4120  ghostGhostedBoundaries(); // We do this again right here in case new boundaries have been added
4121 
4122  // do not assemble system matrix for JFNK solve
4124  _nl->turnOffJacobian();
4125 
4126  {
4127  TIME_SECTION(_eq_init_timer);
4128  _eq.init();
4129  }
4130 
4131  _mesh.meshChanged();
4132  if (_displaced_problem)
4134 
4135  _nl->update();
4136 
4137  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
4138  _assembly[tid]->init(_cm.get());
4139 
4140  _nl->init();
4141 
4142  if (_displaced_problem)
4143  _displaced_problem->init();
4144 
4145  _aux->init();
4146 
4147  _initialized = true;
4148 }
4149 
4150 void
4152 {
4153  TIME_SECTION(_solve_timer);
4154 
4155 #ifdef LIBMESH_HAVE_PETSC
4156  Moose::PetscSupport::petscSetOptions(*this); // Make sure the PETSc options are setup for this app
4157 #endif
4158 
4159  Moose::setSolverDefaults(*this);
4160 
4161  // Setup the output system for printing linear/nonlinear iteration information
4162  initPetscOutput();
4163 
4165 
4166  // reset flag so that linear solver does not use
4167  // the old converged reason "DIVERGED_NANORINF", when
4168  // we throw an exception and stop solve
4170 
4171  if (_solve)
4172  _nl->solve();
4173 
4174  if (_solve)
4175  _nl->update();
4176 
4177  // sync solutions in displaced problem
4178  if (_displaced_problem)
4179  _displaced_problem->syncSolutions();
4180 }
4181 
4182 void
4183 FEProblemBase::setException(const std::string & message)
4184 {
4185  _has_exception = true;
4186  _exception_message = message;
4187 }
4188 
4189 void
4191 {
4193 
4194  // See if any processor had an exception. If it did, get back the
4195  // processor that the exception occurred on.
4196  unsigned int processor_id;
4197 
4198  _communicator.maxloc(_has_exception, processor_id);
4199 
4200  if (_has_exception)
4201  {
4202  _communicator.broadcast(_exception_message, processor_id);
4203 
4204  // Print the message
4205  if (_communicator.rank() == 0)
4206  Moose::err << _exception_message << std::endl;
4207 
4208  // Stop the solve -- this entails setting
4209  // SNESSetFunctionDomainError() or directly inserting NaNs in the
4210  // residual vector to let PETSc >= 3.6 return DIVERGED_NANORINF.
4211  _nl->stopSolve();
4212 
4213  // and close Aux system (we MUST do this here; see #11525)
4214  _aux->solution().close();
4215 
4216  // We've handled this exception, so we no longer have one.
4217  _has_exception = false;
4218 
4219  // Force the next linear convergence check to fail.
4221 
4222  // Repropagate the exception, so it can be caught at a higher level, typically
4223  // this is NonlinearSystem::computeResidual().
4225  }
4226 }
4227 
4228 bool
4230 {
4231  if (_solve)
4232  return _nl->converged();
4233  else
4234  return true;
4235 }
4236 
4237 unsigned int
4239 {
4240  return _nl->nNonlinearIterations();
4241 }
4242 
4243 unsigned int
4245 {
4246  return _nl->nLinearIterations();
4247 }
4248 
4249 Real
4251 {
4252  return _nl->finalNonlinearResidual();
4253 }
4254 
4255 bool
4257 {
4258  return _nl->computingInitialResidual();
4259 }
4260 
4261 void
4263 {
4264  _nl->copySolutionsBackwards();
4265  _aux->copySolutionsBackwards();
4266 }
4267 
4268 void
4270 {
4271  TIME_SECTION(_advance_state_timer);
4272 
4273  _nl->copyOldSolutions();
4274  _aux->copyOldSolutions();
4275 
4276  if (_displaced_problem)
4277  {
4278  _displaced_problem->nlSys().copyOldSolutions();
4279  _displaced_problem->auxSys().copyOldSolutions();
4280  }
4281 
4284 
4286  _material_props.shift(*this);
4287 
4289  _bnd_material_props.shift(*this);
4290 
4293 }
4294 
4295 void
4297 {
4298  TIME_SECTION(_restore_solutions_timer);
4299 
4300  _nl->restoreSolutions();
4301  _aux->restoreSolutions();
4302 
4303  if (_displaced_problem)
4304  _displaced_problem->updateMesh();
4305 }
4306 
4307 void
4309 {
4310  TIME_SECTION(_save_old_solutions_timer);
4311 
4312  _nl->saveOldSolutions();
4313  _aux->saveOldSolutions();
4314 }
4315 
4316 void
4318 {
4319  TIME_SECTION(_restore_old_solutions_timer);
4320 
4321  _nl->restoreOldSolutions();
4322  _aux->restoreOldSolutions();
4323 }
4324 
4325 void
4327 {
4328  TIME_SECTION(_output_step_timer);
4329 
4330  _nl->update();
4331  _aux->update();
4332  if (_displaced_problem)
4333  _displaced_problem->syncSolutions();
4335 }
4336 
4337 void
4339 {
4341 }
4342 
4343 void
4345 {
4347 }
4348 
4349 void
4351 {
4354 }
4355 
4356 void
4358 {
4359  TIME_SECTION(_on_timestep_begin_timer);
4360 
4361  _nl->onTimestepBegin();
4362 }
4363 
4364 void
4366 {
4367 }
4368 
4369 void
4371  const std::string & name,
4372  InputParameters parameters)
4373 {
4374  parameters.set<SubProblem *>("_subproblem") = this;
4375  _aux->addTimeIntegrator(type, name + ":aux", parameters);
4376  _nl->addTimeIntegrator(type, name, parameters);
4377  _has_time_integrator = true;
4378 
4379  // add vectors to store u_dot, u_dotdot, udot_old and u_dotdot_old if requested by the time
4380  // integrator
4381  _aux->addDotVectors();
4382  _nl->addDotVectors();
4383 }
4384 
4385 void
4386 FEProblemBase::addPredictor(const std::string & type,
4387  const std::string & name,
4388  InputParameters parameters)
4389 {
4390  parameters.set<SubProblem *>("_subproblem") = this;
4391  std::shared_ptr<Predictor> predictor = _factory.create<Predictor>(type, name, parameters);
4392  _nl->setPredictor(predictor);
4393 }
4394 
4395 Real
4397 {
4398  TIME_SECTION(_compute_residual_l2_norm_timer);
4399 
4400  computeResidual(*_nl->currentSolution(), _nl->RHS());
4401 
4402  return _nl->RHS().l2_norm();
4403 }
4404 
4405 void
4406 FEProblemBase::computeResidualSys(NonlinearImplicitSystem & /*sys*/,
4407  const NumericVector<Number> & soln,
4408  NumericVector<Number> & residual)
4409 {
4410  TIME_SECTION(_compute_residual_sys_timer);
4411 
4412  computeResidual(soln, residual);
4413 }
4414 
4415 void
4416 FEProblemBase::computeResidual(NonlinearImplicitSystem & sys,
4417  const NumericVector<Number> & soln,
4418  NumericVector<Number> & residual)
4419 {
4420  mooseDeprecated("Please use computeResidualSys");
4421 
4422  computeResidualSys(sys, soln, residual);
4423 }
4424 
4425 void
4426 FEProblemBase::computeResidual(const NumericVector<Number> & soln, NumericVector<Number> & residual)
4427 {
4428  auto & tags = getVectorTags();
4429 
4430  _fe_vector_tags.clear();
4431 
4432  for (auto & tag : tags)
4433  _fe_vector_tags.insert(tag.second);
4434 
4435  computeResidualInternal(soln, residual, _fe_vector_tags);
4436 }
4437 
4438 void
4439 FEProblemBase::computeResidualTag(const NumericVector<Number> & soln,
4440  NumericVector<Number> & residual,
4441  TagID tag)
4442 {
4443  try
4444  {
4445  _nl->setSolution(soln);
4446 
4447  _nl->associateVectorToTag(residual, tag);
4448 
4449  computeResidualTags({tag});
4450 
4451  _nl->disassociateVectorFromTag(residual, tag);
4452  }
4453  catch (MooseException & e)
4454  {
4455  // If a MooseException propagates all the way to here, it means
4456  // that it was thrown from a MOOSE system where we do not
4457  // (currently) properly support the throwing of exceptions, and
4458  // therefore we have no choice but to error out. It may be
4459  // *possible* to handle exceptions from other systems, but in the
4460  // meantime, we don't want to silently swallow any unhandled
4461  // exceptions here.
4462  mooseError("An unhandled MooseException was raised during residual computation. Please "
4463  "contact the MOOSE team for assistance.");
4464  }
4465 }
4466 
4467 void
4468 FEProblemBase::computeResidualInternal(const NumericVector<Number> & soln,
4469  NumericVector<Number> & residual,
4470  const std::set<TagID> & tags)
4471 {
4472  TIME_SECTION(_compute_residual_internal_timer);
4473 
4474  try
4475  {
4476  _nl->setSolution(soln);
4477 
4478  _nl->associateVectorToTag(residual, _nl->residualVectorTag());
4479 
4480  computeResidualTags(tags);
4481 
4482  _nl->disassociateVectorFromTag(residual, _nl->residualVectorTag());
4483  }
4484  catch (MooseException & e)
4485  {
4486  // If a MooseException propagates all the way to here, it means
4487  // that it was thrown from a MOOSE system where we do not
4488  // (currently) properly support the throwing of exceptions, and
4489  // therefore we have no choice but to error out. It may be
4490  // *possible* to handle exceptions from other systems, but in the
4491  // meantime, we don't want to silently swallow any unhandled
4492  // exceptions here.
4493  mooseError("An unhandled MooseException was raised during residual computation. Please "
4494  "contact the MOOSE team for assistance.");
4495  }
4496 }
4497 
4498 void
4499 FEProblemBase::computeResidualType(const NumericVector<Number> & soln,
4500  NumericVector<Number> & residual,
4501  TagID tag)
4502 {
4503  TIME_SECTION(_compute_residual_type_timer);
4504 
4505  try
4506  {
4507  _nl->setSolution(soln);
4508 
4509  _nl->associateVectorToTag(residual, _nl->residualVectorTag());
4510 
4511  computeResidualTags({tag, _nl->residualVectorTag()});
4512 
4513  _nl->disassociateVectorFromTag(residual, _nl->residualVectorTag());
4514  }
4515  catch (MooseException & e)
4516  {
4517  // If a MooseException propagates all the way to here, it means
4518  // that it was thrown from a MOOSE system where we do not
4519  // (currently) properly support the throwing of exceptions, and
4520  // therefore we have no choice but to error out. It may be
4521  // *possible* to handle exceptions from other systems, but in the
4522  // meantime, we don't want to silently swallow any unhandled
4523  // exceptions here.
4524  mooseError("An unhandled MooseException was raised during residual computation. Please "
4525  "contact the MOOSE team for assistance.");
4526  }
4527 }
4528 
4529 void
4531  const NumericVector<Number> & u,
4532  const NumericVector<Number> & udot,
4533  const NumericVector<Number> & udotdot,
4534  NumericVector<Number> & residual)
4535 {
4537 
4538  if (uDotRequested())
4539  _nl->setSolutionUDot(udot);
4540 
4541  if (uDotDotRequested())
4542  _nl->setSolutionUDotDot(udotdot);
4543 
4544  _time = time;
4545  computeResidual(u, residual);
4546 }
4547 
4548 void
4549 FEProblemBase::computeResidualTags(const std::set<TagID> & tags)
4550 {
4551  TIME_SECTION(_compute_residual_tags_timer);
4552 
4553  _nl->zeroVariablesForResidual();
4554  _aux->zeroVariablesForResidual();
4555 
4556  unsigned int n_threads = libMesh::n_threads();
4557 
4559 
4560  // Random interface objects
4561  for (const auto & it : _random_data_objects)
4562  it.second->updateSeeds(EXEC_LINEAR);
4563 
4565 
4567 
4568  for (unsigned int tid = 0; tid < n_threads; tid++)
4569  reinitScalars(tid);
4570 
4572 
4573  _aux->residualSetup();
4574 
4575  if (_displaced_problem)
4576  {
4577  _aux->compute(EXEC_PRE_DISPLACE);
4578  _displaced_problem->updateMesh();
4580  updateMortarMesh();
4581  }
4582 
4583  for (THREAD_ID tid = 0; tid < n_threads; tid++)
4584  {
4587  }
4588 
4589  _nl->computeTimeDerivatives();
4590 
4591  try
4592  {
4593  _aux->compute(EXEC_LINEAR);
4594  }
4595  catch (MooseException & e)
4596  {
4597  _console << "\nA MooseException was raised during Auxiliary variable computation.\n"
4598  << "The next solve will fail, the timestep will be reduced, and we will try again.\n"
4599  << std::endl;
4600 
4601  // We know the next solve is going to fail, so there's no point in
4602  // computing anything else after this. Plus, using incompletely
4603  // computed AuxVariables in subsequent calculations could lead to
4604  // other errors or unhandled exceptions being thrown.
4605  return;
4606  }
4607 
4609 
4611 
4613 
4615 
4617 
4618  _nl->computeResidualTags(tags);
4619 
4621 }
4622 
4623 void
4624 FEProblemBase::computeJacobianSys(NonlinearImplicitSystem & /*sys*/,
4625  const NumericVector<Number> & soln,
4626  SparseMatrix<Number> & jacobian)
4627 {
4628  computeJacobian(soln, jacobian);
4629 }
4630 
4631 void
4632 FEProblemBase::computeJacobianTag(const NumericVector<Number> & soln,
4633  SparseMatrix<Number> & jacobian,
4634  TagID tag)
4635 {
4636  _nl->setSolution(soln);
4637 
4638  _nl->associateMatrixToTag(jacobian, tag);
4639 
4640  computeJacobianTags({tag});
4641 
4642  _nl->disassociateMatrixFromTag(jacobian, tag);
4643 }
4644 
4645 void
4646 FEProblemBase::computeJacobian(const NumericVector<Number> & soln, SparseMatrix<Number> & jacobian)
4647 {
4648  _fe_matrix_tags.clear();
4649 
4650  auto & tags = getMatrixTags();
4651  for (auto & tag : tags)
4652  _fe_matrix_tags.insert(tag.second);
4653 
4654  computeJacobianInternal(soln, jacobian, _fe_matrix_tags);
4655 }
4656 
4657 void
4658 FEProblemBase::computeJacobianInternal(const NumericVector<Number> & soln,
4659  SparseMatrix<Number> & jacobian,
4660  const std::set<TagID> & tags)
4661 {
4662  TIME_SECTION(_compute_jacobian_internal_timer);
4663 
4664  _nl->setSolution(soln);
4665 
4666  _nl->associateMatrixToTag(jacobian, _nl->systemMatrixTag());
4667 
4668  computeJacobianTags(tags);
4669 
4670  _nl->disassociateMatrixFromTag(jacobian, _nl->systemMatrixTag());
4671 }
4672 
4673 void
4674 FEProblemBase::computeJacobianTags(const std::set<TagID> & tags)
4675 {
4676  if (!_has_jacobian || !_const_jacobian)
4677  {
4678  TIME_SECTION(_compute_jacobian_tags_timer);
4679 
4680  for (auto tag : tags)
4681  if (_nl->hasMatrix(tag))
4682  _nl->getMatrix(tag).zero();
4683 
4684  _nl->zeroVariablesForJacobian();
4685  _aux->zeroVariablesForJacobian();
4686 
4687  unsigned int n_threads = libMesh::n_threads();
4688 
4689  // Random interface objects
4690  for (const auto & it : _random_data_objects)
4691  it.second->updateSeeds(EXEC_NONLINEAR);
4692 
4695  if (_displaced_problem)
4696  _displaced_problem->setCurrentlyComputingJacobian(true);
4697 
4700 
4701  for (unsigned int tid = 0; tid < n_threads; tid++)
4702  reinitScalars(tid);
4703 
4705 
4706  _aux->jacobianSetup();
4707 
4708  if (_displaced_problem)
4709  {
4710  _aux->compute(EXEC_PRE_DISPLACE);
4711  _displaced_problem->updateMesh();
4712  }
4713 
4714  for (unsigned int tid = 0; tid < n_threads; tid++)
4715  {
4718  }
4719 
4720  _aux->compute(EXEC_NONLINEAR);
4721 
4723 
4725 
4727 
4729 
4730  _nl->computeJacobianTags(tags);
4731 
4734  if (_displaced_problem)
4735  _displaced_problem->setCurrentlyComputingJacobian(false);
4736  _has_jacobian = true;
4738  }
4739 }
4740 
4741 void
4743  const NumericVector<Number> & u,
4744  const NumericVector<Number> & udot,
4745  const NumericVector<Number> & udotdot,
4746  Real duDotDu_shift,
4747  Real duDotDotDu_shift,
4748  SparseMatrix<Number> & jacobian)
4749 {
4750  if (0)
4751  { // The current interface guarantees that the residual is called before Jacobian, thus udot has
4752  // already been set
4753  if (uDotDotRequested())
4754  _nl->setSolutionUDotDot(udotdot);
4755  if (uDotOldRequested())
4756  _nl->setSolutionUDot(udot);
4757  }
4758  _nl->duDotDu() = duDotDu_shift;
4759  _nl->duDotDotDu() = duDotDotDu_shift;
4760  _time = time;
4761  computeJacobian(u, jacobian);
4762 }
4763 
4764 void
4765 FEProblemBase::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
4766 {
4767  TIME_SECTION(_compute_jacobian_blocks_timer);
4768 
4769  if (_displaced_problem)
4770  {
4771  _aux->compute(EXEC_PRE_DISPLACE);
4772  _displaced_problem->updateMesh();
4773  }
4774 
4775  _aux->compute(EXEC_NONLINEAR);
4776 
4777  _nl->computeJacobianBlocks(blocks);
4778 }
4779 
4780 void
4781 FEProblemBase::computeJacobianBlock(SparseMatrix<Number> & jacobian,
4782  libMesh::System & precond_system,
4783  unsigned int ivar,
4784  unsigned int jvar)
4785 {
4786  JacobianBlock jac_block(precond_system, jacobian, ivar, jvar);
4787  std::vector<JacobianBlock *> blocks = {&jac_block};
4788  computeJacobianBlocks(blocks);
4789 }
4790 
4791 void
4792 FEProblemBase::computeBounds(NonlinearImplicitSystem & /*sys*/,
4793  NumericVector<Number> & lower,
4794  NumericVector<Number> & upper)
4795 {
4796  if (!_nl->hasVector("lower_bound") || !_nl->hasVector("upper_bound"))
4797  return;
4798 
4799  TIME_SECTION(_compute_bounds_timer);
4800 
4801  NumericVector<Number> & _lower = _nl->getVector("lower_bound");
4802  NumericVector<Number> & _upper = _nl->getVector("upper_bound");
4803  _lower.swap(lower);
4804  _upper.swap(upper);
4805  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
4807 
4808  _aux->residualSetup();
4809  _aux->compute(EXEC_LINEAR);
4810  _lower.swap(lower);
4811  _upper.swap(upper);
4812 
4814 }
4815 
4816 void
4817 FEProblemBase::computeNearNullSpace(NonlinearImplicitSystem & /*sys*/,
4818  std::vector<NumericVector<Number> *> & sp)
4819 {
4820  sp.clear();
4821  for (unsigned int i = 0; i < subspaceDim("NearNullSpace"); ++i)
4822  {
4823  std::stringstream postfix;
4824  postfix << "_" << i;
4825  std::string modename = "NearNullSpace" + postfix.str();
4826  sp.push_back(&_nl->getVector(modename));
4827  }
4828 }
4829 
4830 void
4831 FEProblemBase::computeNullSpace(NonlinearImplicitSystem & /*sys*/,
4832  std::vector<NumericVector<Number> *> & sp)
4833 {
4834  sp.clear();
4835  for (unsigned int i = 0; i < subspaceDim("NullSpace"); ++i)
4836  {
4837  std::stringstream postfix;
4838  postfix << "_" << i;
4839  sp.push_back(&_nl->getVector("NullSpace" + postfix.str()));
4840  }
4841 }
4842 
4843 void
4844 FEProblemBase::computeTransposeNullSpace(NonlinearImplicitSystem & /*sys*/,
4845  std::vector<NumericVector<Number> *> & sp)
4846 {
4847  sp.clear();
4848  for (unsigned int i = 0; i < subspaceDim("TransposeNullSpace"); ++i)
4849  {
4850  std::stringstream postfix;
4851  postfix << "_" << i;
4852  sp.push_back(&_nl->getVector("TransposeNullSpace" + postfix.str()));
4853  }
4854 }
4855 
4856 void
4857 FEProblemBase::computePostCheck(NonlinearImplicitSystem & sys,
4858  const NumericVector<Number> & old_soln,
4859  NumericVector<Number> & search_direction,
4860  NumericVector<Number> & new_soln,
4861  bool & changed_search_direction,
4862  bool & changed_new_soln)
4863 {
4864  // This function replaces the old PetscSupport::dampedCheck() function.
4865  //
4866  // 1.) Recreate code in PetscSupport::dampedCheck() for constructing
4867  // ghosted "soln" and "update" vectors.
4868  // 2.) Call FEProblemBase::computeDamping() with these ghost vectors.
4869  // 3.) Recreate the code in PetscSupport::dampedCheck() to actually update
4870  // the solution vector based on the damping, and set the "changed" flags
4871  // appropriately.
4872 
4873  TIME_SECTION(_compute_post_check_timer);
4874 
4875  // MOOSE's FEProblemBase doesn't update the solution during the
4876  // postcheck, but FEProblemBase-derived classes might.
4878  {
4879  // We need ghosted versions of new_soln and search_direction (the
4880  // ones we get from libmesh/PETSc are PARALLEL vectors. To make
4881  // our lives simpler, we use the same ghosting pattern as the
4882  // system's current_local_solution to create new ghosted vectors.
4883 
4884  // Construct zeroed-out clones with the same ghosted dofs as the
4885  // System's current_local_solution.
4886  std::unique_ptr<NumericVector<Number>> ghosted_solution =
4887  sys.current_local_solution->zero_clone(),
4888  ghosted_search_direction =
4889  sys.current_local_solution->zero_clone();
4890 
4891  // Copy values from input vectors into clones with ghosted values.
4892  *ghosted_solution = new_soln;
4893  *ghosted_search_direction = search_direction;
4894 
4895  if (_has_dampers)
4896  {
4897  // Compute the damping coefficient using the ghosted vectors
4898  Real damping = computeDamping(*ghosted_solution, *ghosted_search_direction);
4899 
4900  // If some non-trivial damping was computed, update the new_soln
4901  // vector accordingly.
4902  if (damping < 1.0)
4903  {
4904  new_soln = old_soln;
4905  new_soln.add(-damping, search_direction);
4906  changed_new_soln = true;
4907  }
4908  }
4909 
4910  if (shouldUpdateSolution())
4911  {
4912  // Update the ghosted copy of the new solution, if necessary.
4913  if (changed_new_soln)
4914  *ghosted_solution = new_soln;
4915 
4916  bool updated_solution = updateSolution(new_soln, *ghosted_solution);
4917  if (updated_solution)
4918  changed_new_soln = true;
4919  }
4920  }
4921 
4923  {
4924  _nl->setPreviousNewtonSolution(old_soln);
4925  _aux->setPreviousNewtonSolution();
4926  }
4927 
4928  // MOOSE doesn't change the search_direction
4929  changed_search_direction = false;
4930 }
4931 
4932 Real
4933 FEProblemBase::computeDamping(const NumericVector<Number> & soln,
4934  const NumericVector<Number> & update)
4935 {
4936  // Default to no damping
4937  Real damping = 1.0;
4938 
4939  if (_has_dampers)
4940  {
4941  TIME_SECTION(_compute_damping_timer);
4942 
4943  // Save pointer to the current solution
4944  const NumericVector<Number> * _saved_current_solution = _nl->currentSolution();
4945 
4946  _nl->setSolution(soln);
4947  // For now, do not re-compute auxiliary variables. Doing so allows a wild solution increment
4948  // to get to the material models, which may not be able to cope with drastically different
4949  // values. Once more complete dependency checking is in place, auxiliary variables (and
4950  // material properties) will be computed as needed by dampers.
4951  // _aux.compute();
4952  damping = _nl->computeDamping(soln, update);
4953 
4954  // restore saved solution
4955  _nl->setSolution(*_saved_current_solution);
4956  }
4957 
4958  return damping;
4959 }
4960 
4961 bool
4963 {
4964  return false;
4965 }
4966 
4967 bool
4968 FEProblemBase::updateSolution(NumericVector<Number> & /*vec_solution*/,
4969  NumericVector<Number> & /*ghosted_solution*/)
4970 {
4971  return false;
4972 }
4973 
4974 void
4975 FEProblemBase::predictorCleanup(NumericVector<Number> & /*ghosted_solution*/)
4976 {
4977 }
4978 
4979 void
4980 FEProblemBase::addDisplacedProblem(std::shared_ptr<DisplacedProblem> displaced_problem)
4981 {
4982  _displaced_mesh = &displaced_problem->mesh();
4983  _displaced_problem = displaced_problem;
4984 }
4985 
4986 void
4988 {
4989  TIME_SECTION(_update_geometric_search_timer);
4990 
4992 
4993  if (_displaced_problem)
4994  _displaced_problem->updateGeomSearch(type);
4995 }
4996 
4997 void
4999 {
5000  _mortar_data.update();
5001 }
5002 
5003 void
5005  const std::pair<BoundaryID, BoundaryID> & master_slave_boundary_pair,
5006  const std::pair<SubdomainID, SubdomainID> & master_slave_subdomain_pair,
5007  bool on_displaced,
5008  bool periodic)
5009 {
5010  if (on_displaced)
5011  return _mortar_data.createMortarInterface(master_slave_boundary_pair,
5012  master_slave_subdomain_pair,
5014  on_displaced,
5015  periodic);
5016  else
5018  master_slave_boundary_pair, master_slave_subdomain_pair, *this, on_displaced, periodic);
5019 }
5020 
5023  const std::pair<BoundaryID, BoundaryID> & master_slave_boundary_pair,
5024  const std::pair<SubdomainID, SubdomainID> & master_slave_subdomain_pair,
5025  bool on_displaced) const
5026 {
5028  master_slave_boundary_pair, master_slave_subdomain_pair, on_displaced);
5029 }
5030 
5031 const std::unordered_map<std::pair<BoundaryID, BoundaryID>, AutomaticMortarGeneration> &
5032 FEProblemBase::getMortarInterfaces(bool on_displaced) const
5033 {
5034  return _mortar_data.getMortarInterfaces(on_displaced);
5035 }
5036 
5037 void
5039 {
5040  if (_displaced_problem) // Only need to do this if things are moving...
5041  {
5043 
5044  switch (_mesh.getPatchUpdateStrategy())
5045  {
5046  case Moose::Never:
5047  break;
5048  case Moose::Iteration:
5049  // Update the list of ghosted elements at the start of the time step
5052 
5053  _displaced_problem->geomSearchData().updateGhostedElems();
5055 
5056  // The commands below ensure that the sparsity of the Jacobian matrix is
5057  // augmented at the start of the time step using neighbor nodes from the end
5058  // of the previous time step.
5059 
5061 
5062  // This is needed to reinitialize PETSc output
5063  initPetscOutput();
5064 
5065  break;
5066 
5067  case Moose::Auto:
5068  {
5069  Real max = _displaced_problem->geomSearchData().maxPatchPercentage();
5070  _communicator.max(max);
5071 
5072  // If we haven't moved very far through the patch
5073  if (max < 0.4)
5074  break;
5075  }
5076  libmesh_fallthrough();
5077 
5078  // Let this fall through if things do need to be updated...
5079  case Moose::Always:
5080  // Flush output here to see the message before the reinitialization, which could take a
5081  // while
5082  _console << "\n\nUpdating geometric search patches\n" << std::endl;
5083 
5086 
5087  _displaced_problem->geomSearchData().clearNearestNodeLocators();
5089 
5091 
5092  // This is needed to reinitialize PETSc output
5093  initPetscOutput();
5094  }
5095  }
5096 }
5097 
5098 #ifdef LIBMESH_ENABLE_AMR
5099 void
5101 {
5102  unsigned int n = adaptivity().getInitialSteps();
5103  _cycles_completed = 0;
5104  if (n)
5105  {
5106  TIME_SECTION(_initial_adapt_mesh_timer);
5107 
5108  for (unsigned int i = 0; i < n; i++)
5109  {
5110  _console << "Initial adaptivity step " << i + 1 << " of " << n << std::endl;
5112  computeMarkers();
5113 
5115  {
5116  meshChanged();
5117 
5118  // reproject the initial condition
5119  projectSolution();
5120 
5122  }
5123  else
5124  {
5125  _console << "Mesh unchanged, skipping remaining steps..." << std::endl;
5126  return;
5127  }
5128  }
5129  }
5130 }
5131 
5132 bool
5134 {
5135  // reset cycle counter
5136  _cycles_completed = 0;
5137 
5139  return false;
5140 
5141  TIME_SECTION(_adapt_mesh_timer);
5142 
5143  unsigned int cycles_per_step = _adaptivity.getCyclesPerStep();
5144 
5145  bool mesh_changed = false;
5146 
5147  for (unsigned int i = 0; i < cycles_per_step; ++i)
5148  {
5149  _console << "Adaptivity step " << i + 1 << " of " << cycles_per_step << '\n';
5150 
5151  // Markers were already computed once by Executioner
5152  if (_adaptivity.getRecomputeMarkersFlag() && i > 0)
5153  computeMarkers();
5154 
5155  if (_adaptivity.adaptMesh())
5156  {
5157  mesh_changed = true;
5158 
5159  meshChangedHelper(true); // This may be an intermediate change
5161  }
5162  else
5163  {
5164  _console << "Mesh unchanged, skipping remaining steps..." << std::endl;
5165  break;
5166  }
5167 
5168  // Show adaptivity progress
5169  _console << std::flush;
5170  }
5171 
5172  // We're done with all intermediate changes; now get systems ready
5173  // for real if necessary.
5174  if (mesh_changed)
5175  _eq.reinit_systems();
5176 
5177  return mesh_changed;
5178 }
5179 #endif // LIBMESH_ENABLE_AMR
5180 
5181 void
5182 FEProblemBase::initXFEM(std::shared_ptr<XFEMInterface> xfem)
5183 {
5184  _xfem = xfem;
5185  _xfem->setMesh(&_mesh);
5186  if (_displaced_mesh)
5187  _xfem->setDisplacedMesh(_displaced_mesh);
5188  _xfem->setMaterialData(&_material_data);
5189  _xfem->setBoundaryMaterialData(&_bnd_material_data);
5190 
5191  unsigned int n_threads = libMesh::n_threads();
5192  for (unsigned int i = 0; i < n_threads; ++i)
5193  {
5194  _assembly[i]->setXFEM(_xfem);
5195  if (_displaced_problem)
5196  _displaced_problem->assembly(i).setXFEM(_xfem);
5197  }
5198 }
5199 
5200 bool
5202 {
5203  TIME_SECTION(_update_mesh_xfem_timer);
5204 
5205  bool updated = false;
5206  if (haveXFEM())
5207  {
5208  if (_xfem->updateHeal())
5209  meshChanged();
5210 
5211  updated = _xfem->update(_time, *_nl, *_aux);
5212  if (updated)
5213  {
5214  meshChanged();
5215  _xfem->initSolution(*_nl, *_aux);
5216  restoreSolutions();
5217  }
5218  }
5219  return updated;
5220 }
5221 
5222 void
5224 {
5225  TIME_SECTION(_mesh_changed_timer);
5226 
5227  this->meshChangedHelper();
5228 }
5229 
5230 void
5231 FEProblemBase::meshChangedHelper(bool intermediate_change)
5232 {
5233  TIME_SECTION(_mesh_changed_helper_timer);
5234 
5237  _mesh.cacheChangedLists(); // Currently only used with adaptivity and stateful material
5238  // properties
5239 
5240  // Clear these out because they corresponded to the old mesh
5241  _ghosted_elems.clear();
5242 
5244 
5245  // The mesh changed. We notify the MooseMesh first, because
5246  // callbacks (e.g. for sparsity calculations) triggered by the
5247  // EquationSystems reinit may require up-to-date MooseMesh caches.
5248  _mesh.meshChanged();
5249 
5250  // If we're just going to alter the mesh again, all we need to
5251  // handle here is AMR and projections, not full system reinit
5252  if (intermediate_change)
5253  _eq.reinit_solutions();
5254  else
5255  _eq.reinit();
5256 
5257  // Updating MooseMesh first breaks other adaptivity code, unless we
5258  // then *again* update the MooseMesh caches. E.g. the definition of
5259  // "active" and "local" may have been *changed* by refinement and
5260  // repartitioning done in EquationSystems::reinit().
5261  _mesh.meshChanged();
5262 
5263  // Since the Mesh changed, update the PointLocator object used by DiracKernels.
5265 
5266  // Need to redo ghosting
5268 
5269  if (_displaced_problem)
5270  {
5271  _displaced_problem->meshChanged();
5273  }
5274 
5276 
5278 
5280 
5281  // We need to create new storage for the new elements and copy stateful properties from the old
5282  // elements.
5285  {
5286  {
5287  ProjectMaterialProperties pmp(true,
5288  *this,
5289  *_nl,
5294  _assembly);
5295  Threads::parallel_reduce(*_mesh.refinedElementRange(), pmp);
5296  }
5297 
5298  {
5299  ProjectMaterialProperties pmp(false,
5300  *this,
5301  *_nl,
5306  _assembly);
5307  Threads::parallel_reduce(*_mesh.coarsenedElementRange(), pmp);
5308  }
5309  }
5310 
5313 
5314  _has_jacobian = false; // we have to recompute jacobian when mesh changed
5315 
5316  for (const auto & mci : _notify_when_mesh_changes)
5317  mci->meshChanged();
5318 }
5319 
5320 void
5322 {
5323  _notify_when_mesh_changes.push_back(mci);
5324 }
5325 
5326 void
5328 {
5329  TIME_SECTION(_check_problem_integrity_timer);
5330 
5331  // Check for unsatisfied actions
5332  const std::set<SubdomainID> & mesh_subdomains = _mesh.meshSubdomains();
5333 
5334  // Check kernel coverage of subdomains (blocks) in the mesh
5336  _nl->checkKernelCoverage(mesh_subdomains);
5337 
5338  // Check materials
5339  {
5340 #ifdef LIBMESH_ENABLE_AMR
5341  if (_adaptivity.isOn() &&
5344  {
5345  _console << "Using EXPERIMENTAL Stateful Material Property projection with Adaptivity!\n";
5346 
5347  if (n_processors() > 1)
5348  {
5349  if (_mesh.uniformRefineLevel() > 0 && _mesh.getMesh().skip_partitioning() == false)
5350  mooseError("This simulation is using uniform refinement on the mesh, with stateful "
5351  "properties and adaptivity. "
5352  "You must skip partitioning to run this case:\nMesh/skip_partitioning=true");
5353 
5354  _console << "\nWarning! Mesh re-partitioning is disabled while using stateful material "
5355  "properties! This can lead to large load imbalances and degraded "
5356  "performance!!\n\n";
5357  _mesh.getMesh().skip_partitioning(true);
5358 
5359  _mesh.errorIfDistributedMesh("StatefulMaterials + Adaptivity");
5360 
5361  if (_displaced_problem)
5362  _displaced_problem->mesh().getMesh().skip_partitioning(true);
5363  }
5364  }
5365 #endif
5366 
5367  std::set<SubdomainID> local_mesh_subs(mesh_subdomains);
5368 
5370  {
5375  bool check_material_coverage = false;
5376  std::set<SubdomainID> ids = _all_materials.getActiveBlocks();
5377  for (const auto & id : ids)
5378  {
5379  local_mesh_subs.erase(id);
5380  check_material_coverage = true;
5381  }
5382 
5383  // also exclude mortar spaces from the material check
5384  auto && mortar_subdomain_ids = _mortar_data.getMortarSubdomainIDs();
5385  for (auto subdomain_id : mortar_subdomain_ids)
5386  local_mesh_subs.erase(subdomain_id);
5387 
5388  // Check Material Coverage
5389  if (check_material_coverage && !local_mesh_subs.empty())
5390  {
5391  std::stringstream extra_subdomain_ids;
5393  std::copy(local_mesh_subs.begin(),
5394  local_mesh_subs.end(),
5395  std::ostream_iterator<unsigned int>(extra_subdomain_ids, " "));
5396 
5397  mooseError("The following blocks from your input mesh do not contain an active material: " +
5398  extra_subdomain_ids.str() +
5399  "\nWhen ANY mesh block contains a Material object, "
5400  "all blocks must contain a Material object.\n");
5401  }
5402  }
5403 
5404  // Check material properties on blocks and boundaries
5407 
5408  // Check that material properties exist when requested by other properties on a given block
5409  const auto & materials = _all_materials.getActiveObjects();
5410  for (const auto & material : materials)
5411  material->checkStatefulSanity();
5412 
5413  // auto mats_to_check = _materials.getActiveBlockObjects();
5414  // const auto & discrete_materials = _discrete_materials.getActiveBlockObjects();
5415  // for (const auto & map_it : discrete_materials)
5416  // for (const auto & container_element : map_it.second)
5417  // mats_to_check[map_it.first].push_back(container_element);
5419  }
5420 
5421  // Check UserObjects and Postprocessors
5422  checkUserObjects();
5423 
5424  // Verify that we don't have any Element type/Coordinate Type conflicts
5426 
5427  // If using displacements, verify that the order of the displacement
5428  // variables matches the order of the elements in the displaced
5429  // mesh.
5431 }
5432 
5433 void
5435 {
5436  if (_displaced_problem)
5437  {
5438  bool mesh_has_second_order_elements = false;
5439  for (const auto & elem : as_range(_displaced_mesh->activeLocalElementsBegin(),
5441  {
5442  if (elem->default_order() == SECOND)
5443  {
5444  mesh_has_second_order_elements = true;
5445  break;
5446  }
5447  }
5448 
5449  // We checked our local elements, so take the max over all processors.
5450  _displaced_mesh->comm().max(mesh_has_second_order_elements);
5451 
5452  // If the Mesh has second order elements, make sure the
5453  // displacement variables are second-order.
5454  if (mesh_has_second_order_elements)
5455  {
5456  const std::vector<std::string> & displacement_variables =
5457  _displaced_problem->getDisplacementVarNames();
5458 
5459  for (const auto & var_name : displacement_variables)
5460  {
5461  MooseVariableFEBase & mv =
5462  _displaced_problem->getVariable(/*tid=*/0,
5463  var_name,
5466  if (mv.order() != SECOND)
5467  mooseError("Error: mesh has SECOND order elements, so all displacement variables must be "
5468  "SECOND order.");
5469  }
5470  }
5471  }
5472 }
5473 
5474 void
5476 {
5477  // Check user_objects block coverage
5478  std::set<SubdomainID> mesh_subdomains = _mesh.meshSubdomains();
5479  std::set<SubdomainID> user_objects_blocks;
5480 
5481  // gather names of all user_objects that were defined in the input file
5482  // and the blocks that they are defined on
5483  std::set<std::string> names;
5484 
5485  std::vector<UserObject *> objects;
5487 
5488  for (const auto & obj : objects)
5489  names.insert(obj->name());
5490 
5491  // See if all referenced blocks are covered
5492  mesh_subdomains.insert(Moose::ANY_BLOCK_ID);
5493  std::set<SubdomainID> difference;
5494  std::set_difference(user_objects_blocks.begin(),
5495  user_objects_blocks.end(),
5496  mesh_subdomains.begin(),
5497  mesh_subdomains.end(),
5498  std::inserter(difference, difference.end()));
5499 
5500  if (!difference.empty())
5501  {
5502  std::ostringstream oss;
5503  oss << "One or more UserObjects is referencing a nonexistent block:\n";
5504  for (const auto & id : difference)
5505  oss << id << "\n";
5506  mooseError(oss.str());
5507  }
5508 
5509  // check that all requested UserObjects were defined in the input file
5510  for (const auto & it : _pps_data.values())
5511  {
5512  if (names.find(it.first) == names.end())
5513  mooseError("Postprocessor '" + it.first + "' requested but not specified in the input file.");
5514  }
5515 }
5516 
5517 void
5519  const std::map<SubdomainID, std::vector<std::shared_ptr<Material>>> & materials_map)
5520 {
5521  auto & prop_names = _material_props.statefulPropNames();
5522 
5523  for (const auto & it : materials_map)
5524  {
5526  std::set<std::string> block_depend_props, block_supplied_props;
5527 
5528  for (const auto & mat1 : it.second)
5529  {
5530  const std::set<std::string> & depend_props = mat1->getRequestedItems();
5531  block_depend_props.insert(depend_props.begin(), depend_props.end());
5532 
5533  auto & alldeps = mat1->getMatPropDependencies(); // includes requested stateful props
5534  for (auto & dep : alldeps)
5535  {
5536  if (prop_names.count(dep) > 0)
5537  block_depend_props.insert(prop_names.at(dep));