Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // MOOSE includes
11 :
12 : #include "AuxiliarySystem.h"
13 : #include "FEProblem.h"
14 : #include "MooseApp.h"
15 : #include "MooseMesh.h"
16 : #include "NonlinearSystem.h"
17 : #include "Problem.h"
18 : #include "ResetDisplacedMeshThread.h"
19 : #include "SubProblem.h"
20 : #include "AllNodesSendListThread.h"
21 : #include "Assembly.h"
22 : #include "DisplacedProblem.h"
23 : #include "libmesh/numeric_vector.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/transient_system.h"
27 : #include "libmesh/explicit_system.h"
28 :
29 : using namespace libMesh;
30 :
31 : registerMooseObject("MooseApp", DisplacedProblem);
32 :
33 : InputParameters
34 7105 : DisplacedProblem::validParams()
35 : {
36 7105 : InputParameters params = SubProblem::validParams();
37 14210 : params.addClassDescription(
38 : "A Problem object for providing access to the displaced finite element "
39 : "mesh and associated variables.");
40 7105 : params.addPrivateParam<MooseMesh *>("mesh");
41 21315 : params.addPrivateParam<std::vector<std::string>>("displacements", {});
42 7105 : return params;
43 0 : }
44 :
45 2022 : DisplacedProblem::DisplacedProblem(const InputParameters & parameters)
46 : : SubProblem(parameters),
47 2022 : _mproblem(parameters.have_parameter<FEProblemBase *>("_fe_problem_base")
48 6066 : ? *getParam<FEProblemBase *>("_fe_problem_base")
49 2022 : : *getParam<FEProblem *>("_fe_problem")),
50 4044 : _mesh(*getParam<MooseMesh *>("mesh")),
51 2022 : _eq(_mesh),
52 2022 : _ref_mesh(_mproblem.mesh()),
53 4044 : _displacements(getParam<std::vector<std::string>>("displacements")),
54 4044 : _geometric_search_data(*this, _mesh)
55 :
56 : {
57 : // Disable refinement/coarsening in EquationSystems::reinit because we already do this ourselves
58 2022 : _eq.disable_refine_in_reinit();
59 :
60 : // TODO: Move newAssemblyArray further up to SubProblem so that we can use it here
61 2022 : unsigned int n_threads = libMesh::n_threads();
62 :
63 2022 : _assembly.resize(n_threads);
64 4044 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
65 : {
66 2022 : _displaced_solver_systems.emplace_back(std::make_unique<DisplacedSystem>(
67 : *this,
68 : _mproblem,
69 2022 : _mproblem.getNonlinearSystemBase(nl_sys_num),
70 4044 : "displaced_" + _mproblem.getNonlinearSystemBase(nl_sys_num).name() + "_" +
71 2022 : std::to_string(nl_sys_num),
72 2022 : Moose::VAR_SOLVER));
73 2022 : auto & displaced_nl = _displaced_solver_systems.back();
74 :
75 4244 : for (unsigned int i = 0; i < n_threads; ++i)
76 2222 : _assembly[i].emplace_back(std::make_unique<Assembly>(*displaced_nl, i));
77 : }
78 :
79 2022 : _nl_solution.resize(_displaced_solver_systems.size(), nullptr);
80 :
81 : _displaced_aux =
82 4044 : std::make_unique<DisplacedSystem>(*this,
83 : _mproblem,
84 2022 : _mproblem.getAuxiliarySystem(),
85 4044 : "displaced_" + _mproblem.getAuxiliarySystem().name(),
86 4044 : Moose::VAR_AUXILIARY);
87 :
88 : // // Generally speaking, the mesh is prepared for use, and consequently remote elements are deleted
89 : // // well before our Problem(s) are constructed. Historically, in MooseMesh we have a bunch of
90 : // // needs_prepare type flags that make it so we never call prepare_for_use (and consequently
91 : // // delete_remote_elements) again. So the below line, historically, has had no impact. HOWEVER:
92 : // // I've added some code in SetupMeshCompleteAction for deleting remote elements post
93 : // // EquationSystems::init. If I execute that code without default ghosting, then I get > 40 MOOSE
94 : // // test failures, so we clearly have some simulations that are not yet covered properly by
95 : // // relationship managers. Until that is resolved, I am going to retain default geometric ghosting
96 : // if (!_default_ghosting)
97 : // _mesh.getMesh().remove_ghosting_functor(_mesh.getMesh().default_ghosting());
98 :
99 2022 : automaticScaling(_mproblem.automaticScaling());
100 :
101 2022 : _mesh.setCoordData(_ref_mesh);
102 2022 : }
103 :
104 1992 : DisplacedProblem::~DisplacedProblem() = default;
105 :
106 : bool
107 92663617 : DisplacedProblem::isTransient() const
108 : {
109 92663617 : return _mproblem.isTransient();
110 : }
111 :
112 : std::set<dof_id_type> &
113 8658 : DisplacedProblem::ghostedElems()
114 : {
115 8658 : return _mproblem.ghostedElems();
116 : }
117 :
118 : void
119 2022 : DisplacedProblem::createQRules(QuadratureType type,
120 : Order order,
121 : Order volume_order,
122 : Order face_order,
123 : SubdomainID block,
124 : const bool allow_negative_qweights)
125 : {
126 4244 : for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
127 4444 : for (const auto sys_num : index_range(_assembly[tid]))
128 2222 : _assembly[tid][sys_num]->createQRules(
129 : type, order, volume_order, face_order, block, allow_negative_qweights);
130 2022 : }
131 :
132 : void
133 0 : DisplacedProblem::bumpVolumeQRuleOrder(Order order, SubdomainID block)
134 : {
135 0 : for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
136 0 : for (const auto nl_sys_num : index_range(_assembly[tid]))
137 0 : _assembly[tid][nl_sys_num]->bumpVolumeQRuleOrder(order, block);
138 0 : }
139 :
140 : void
141 0 : DisplacedProblem::bumpAllQRuleOrder(Order order, SubdomainID block)
142 : {
143 0 : for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
144 0 : for (const auto nl_sys_num : index_range(_assembly[tid]))
145 0 : _assembly[tid][nl_sys_num]->bumpAllQRuleOrder(order, block);
146 0 : }
147 :
148 : void
149 2022 : DisplacedProblem::init()
150 : {
151 4243 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
152 : {
153 4442 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
154 2221 : _assembly[tid][nl_sys_num]->init(_mproblem.couplingMatrix(nl_sys_num));
155 :
156 4442 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
157 : {
158 2221 : std::vector<std::pair<unsigned int, unsigned short>> disp_numbers_and_directions;
159 7027 : for (const auto direction : index_range(_displacements))
160 : {
161 4806 : const auto & disp_string = _displacements[direction];
162 4806 : const auto & disp_variable = getVariable(tid, disp_string);
163 4806 : if (disp_variable.sys().number() == nl_sys_num)
164 2157 : disp_numbers_and_directions.push_back(
165 4314 : std::make_pair(disp_variable.number(), cast_int<unsigned short>(direction)));
166 : }
167 2221 : _assembly[tid][nl_sys_num]->assignDisplacements(std::move(disp_numbers_and_directions));
168 2221 : }
169 : }
170 :
171 4044 : for (auto & nl : _displaced_solver_systems)
172 : {
173 2022 : nl->dofMap().attach_extra_send_list_function(&extraSendList, nl.get());
174 2022 : nl->preInit();
175 : }
176 :
177 2022 : _displaced_aux->dofMap().attach_extra_send_list_function(&extraSendList, _displaced_aux.get());
178 2022 : _displaced_aux->preInit();
179 :
180 : {
181 10110 : TIME_SECTION("eq::init", 2, "Initializing Displaced Equation System");
182 2022 : _eq.init();
183 2022 : }
184 :
185 4044 : for (auto & nl : _displaced_solver_systems)
186 2022 : nl->postInit();
187 2022 : _displaced_aux->postInit();
188 :
189 2022 : _mesh.meshChanged();
190 :
191 2022 : if (haveFV())
192 12 : _mesh.setupFiniteVolumeMeshData();
193 2022 : }
194 :
195 : void
196 124 : DisplacedProblem::initAdaptivity()
197 : {
198 124 : }
199 :
200 : void
201 1728 : DisplacedProblem::addTimeIntegrator()
202 : {
203 3456 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
204 1728 : _displaced_solver_systems[nl_sys_num]->copyTimeIntegrators(
205 1728 : _mproblem.getNonlinearSystemBase(nl_sys_num));
206 1728 : _displaced_aux->copyTimeIntegrators(_mproblem.getAuxiliarySystem());
207 1728 : }
208 :
209 : void
210 10 : DisplacedProblem::saveOldSolutions()
211 : {
212 20 : for (auto & displaced_nl : _displaced_solver_systems)
213 10 : displaced_nl->saveOldSolutions();
214 10 : _displaced_aux->saveOldSolutions();
215 10 : }
216 :
217 : void
218 10 : DisplacedProblem::restoreOldSolutions()
219 : {
220 20 : for (auto & displaced_nl : _displaced_solver_systems)
221 10 : displaced_nl->restoreOldSolutions();
222 10 : _displaced_aux->restoreOldSolutions();
223 10 : }
224 :
225 : void
226 488396 : DisplacedProblem::syncAuxSolution(const NumericVector<Number> & aux_soln)
227 : {
228 488396 : (*_displaced_aux->sys().solution) = aux_soln;
229 488396 : _displaced_aux->update();
230 488396 : }
231 :
232 : void
233 308790 : DisplacedProblem::syncSolutions()
234 : {
235 1543950 : TIME_SECTION("syncSolutions", 5, "Syncing Displaced Solutions");
236 :
237 617580 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
238 : {
239 308790 : auto & displaced_nl = _displaced_solver_systems[nl_sys_num];
240 : mooseAssert(nl_sys_num == displaced_nl->number(),
241 : "We should have designed things such that the nl system numbers make their system "
242 : "numbering in the EquationSystems object");
243 308790 : (*displaced_nl->sys().solution) =
244 617580 : *_mproblem.getNonlinearSystemBase(displaced_nl->number()).currentSolution();
245 308790 : displaced_nl->update();
246 : }
247 308790 : syncAuxSolution(*_mproblem.getAuxiliarySystem().currentSolution());
248 308790 : }
249 :
250 : void
251 0 : DisplacedProblem::syncSolutions(
252 : const std::map<unsigned int, const NumericVector<Number> *> & nl_solns,
253 : const NumericVector<Number> & aux_soln)
254 : {
255 0 : TIME_SECTION("syncSolutions", 5, "Syncing Displaced Solutions");
256 :
257 0 : for (const auto [nl_sys_num, nl_soln] : nl_solns)
258 : {
259 0 : (*_displaced_solver_systems[nl_sys_num]->sys().solution) = *nl_soln;
260 0 : _displaced_solver_systems[nl_sys_num]->update();
261 : }
262 0 : syncAuxSolution(aux_soln);
263 0 : }
264 :
265 : void
266 148216 : DisplacedProblem::updateMesh(bool mesh_changing)
267 : {
268 741080 : TIME_SECTION("updateMesh", 3, "Updating Displaced Mesh");
269 :
270 : // If the mesh is changing, we are probably performing adaptivity. In that case, we do *not* want
271 : // to use the undisplaced mesh solution because it may be out-of-sync, whereas our displaced mesh
272 : // solution should be in the correct state after getting restricted/prolonged in
273 : // EquationSystems::reinit (must have been called before this method)
274 148216 : if (!mesh_changing)
275 147647 : syncSolutions();
276 :
277 296432 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
278 148216 : _nl_solution[nl_sys_num] = _displaced_solver_systems[nl_sys_num]->sys().solution.get();
279 148216 : _aux_solution = _displaced_aux->sys().solution.get();
280 :
281 : // If the displaced mesh has been serialized to one processor (as
282 : // may have occurred if it was used for Exodus output), then we need
283 : // the reference mesh to be also. For that matter, did anyone
284 : // somehow serialize the whole mesh? Hopefully not but let's avoid
285 : // causing errors if so.
286 148216 : if (_mesh.getMesh().is_serial() && !this->refMesh().getMesh().is_serial())
287 0 : this->refMesh().getMesh().allgather();
288 :
289 148216 : if (_mesh.getMesh().is_serial_on_zero() && !this->refMesh().getMesh().is_serial_on_zero())
290 0 : this->refMesh().getMesh().gather_to_zero();
291 :
292 148216 : UpdateDisplacedMeshThread udmt(_mproblem, *this);
293 :
294 : // We displace all nodes, not just semilocal nodes, because
295 : // parallel-inconsistent mesh geometry makes libMesh cry.
296 296432 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
297 296432 : _mesh.getMesh().nodes_end(),
298 148216 : /*grainsize=*/1);
299 :
300 148216 : Threads::parallel_reduce(node_range, udmt);
301 : // Displacement of the mesh has invalidated the point locator data (e.g. bounding boxes)
302 148216 : _mesh.getMesh().clear_point_locator();
303 :
304 : // The mesh has changed. Face information normals, areas, etc. must be re-calculated
305 148216 : if (haveFV())
306 12 : _mesh.setupFiniteVolumeMeshData();
307 :
308 : // Update the geometric searches that depend on the displaced mesh. This call can end up running
309 : // NearestNodeThread::operator() which has a throw inside of it. We need to catch it and make sure
310 : // it's propagated to all processes before updating the point locator because the latter requires
311 : // communication
312 : try
313 : {
314 : // We may need to re-run geometric operations like SecondaryNeighborhoodTread if, for instance,
315 : // we have performed mesh adaptivity
316 148216 : if (mesh_changing)
317 569 : _geometric_search_data.reinit();
318 : else
319 147647 : _geometric_search_data.update();
320 : }
321 0 : catch (MooseException & e)
322 : {
323 0 : _mproblem.setException(e.what());
324 0 : }
325 :
326 148213 : if (udmt.hasDisplacement())
327 63048 : _mproblem.meshDisplaced();
328 :
329 : // The below call will throw an exception on all processes if any of our processes had an
330 : // exception above. This exception will be caught higher up the call stack and the error message
331 : // will be printed there
332 148213 : _mproblem.checkExceptionAndStopSolve(/*print_message=*/false);
333 :
334 : // Since the Mesh changed, update the PointLocator object used by DiracKernels.
335 148213 : _dirac_kernel_info.updatePointLocator(_mesh);
336 148213 : }
337 :
338 : void
339 0 : DisplacedProblem::updateMesh(const std::map<unsigned int, const NumericVector<Number> *> & nl_solns,
340 : const NumericVector<Number> & aux_soln)
341 : {
342 0 : TIME_SECTION("updateMesh", 3, "Updating Displaced Mesh");
343 :
344 0 : syncSolutions(nl_solns, aux_soln);
345 :
346 0 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
347 0 : _nl_solution[nl_sys_num] = _displaced_solver_systems[nl_sys_num]->sys().solution.get();
348 0 : _aux_solution = _displaced_aux->sys().solution.get();
349 :
350 0 : UpdateDisplacedMeshThread udmt(_mproblem, *this);
351 :
352 : // We displace all nodes, not just semilocal nodes, because
353 : // parallel-inconsistent mesh geometry makes libMesh cry.
354 0 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
355 0 : _mesh.getMesh().nodes_end(),
356 0 : /*grainsize=*/1);
357 :
358 0 : Threads::parallel_reduce(node_range, udmt);
359 :
360 : // Update the geometric searches that depend on the displaced mesh. This call can end up running
361 : // NearestNodeThread::operator() which has a throw inside of it. We need to catch it and make sure
362 : // it's propagated to all processes before updating the point locator because the latter requires
363 : // communication
364 : try
365 : {
366 0 : _geometric_search_data.update();
367 : }
368 0 : catch (MooseException & e)
369 : {
370 0 : _mproblem.setException(e.what());
371 0 : }
372 :
373 0 : if (udmt.hasDisplacement())
374 0 : _mproblem.meshDisplaced();
375 :
376 : // The below call will throw an exception on all processes if any of our processes had an
377 : // exception above. This exception will be caught higher up the call stack and the error message
378 : // will be printed there
379 0 : _mproblem.checkExceptionAndStopSolve(/*print_message=*/false);
380 :
381 : // Since the Mesh changed, update the PointLocator object used by DiracKernels.
382 0 : _dirac_kernel_info.updatePointLocator(_mesh);
383 0 : }
384 :
385 : TagID
386 0 : DisplacedProblem::addVectorTag(const TagName & tag_name,
387 : const Moose::VectorTagType type /* = Moose::VECTOR_TAG_RESIDUAL */)
388 : {
389 0 : return _mproblem.addVectorTag(tag_name, type);
390 : }
391 :
392 : const VectorTag &
393 5396986 : DisplacedProblem::getVectorTag(const TagID tag_id) const
394 : {
395 5396986 : return _mproblem.getVectorTag(tag_id);
396 : }
397 :
398 : TagID
399 106475 : DisplacedProblem::getVectorTagID(const TagName & tag_name) const
400 : {
401 106475 : return _mproblem.getVectorTagID(tag_name);
402 : }
403 :
404 : TagName
405 0 : DisplacedProblem::vectorTagName(const TagID tag_id) const
406 : {
407 0 : return _mproblem.vectorTagName(tag_id);
408 : }
409 :
410 : bool
411 213 : DisplacedProblem::vectorTagExists(const TagID tag_id) const
412 : {
413 213 : return _mproblem.vectorTagExists(tag_id);
414 : }
415 :
416 : bool
417 0 : DisplacedProblem::vectorTagExists(const TagName & tag_name) const
418 : {
419 0 : return _mproblem.vectorTagExists(tag_name);
420 : }
421 :
422 : unsigned int
423 44613 : DisplacedProblem::numVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
424 : {
425 44613 : return _mproblem.numVectorTags(type);
426 : }
427 :
428 : const std::vector<VectorTag> &
429 2222 : DisplacedProblem::getVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
430 : {
431 2222 : return _mproblem.getVectorTags(type);
432 : }
433 :
434 : Moose::VectorTagType
435 95089710 : DisplacedProblem::vectorTagType(const TagID tag_id) const
436 : {
437 95089710 : return _mproblem.vectorTagType(tag_id);
438 : }
439 :
440 : TagID
441 0 : DisplacedProblem::addMatrixTag(TagName tag_name)
442 : {
443 0 : return _mproblem.addMatrixTag(tag_name);
444 : }
445 :
446 : TagID
447 1512 : DisplacedProblem::getMatrixTagID(const TagName & tag_name) const
448 : {
449 1512 : return _mproblem.getMatrixTagID(tag_name);
450 : }
451 :
452 : TagName
453 0 : DisplacedProblem::matrixTagName(TagID tag)
454 : {
455 0 : return _mproblem.matrixTagName(tag);
456 : }
457 :
458 : bool
459 0 : DisplacedProblem::matrixTagExists(const TagName & tag_name) const
460 : {
461 0 : return _mproblem.matrixTagExists(tag_name);
462 : }
463 :
464 : bool
465 349 : DisplacedProblem::matrixTagExists(TagID tag_id) const
466 : {
467 349 : return _mproblem.matrixTagExists(tag_id);
468 : }
469 :
470 : unsigned int
471 2221 : DisplacedProblem::numMatrixTags() const
472 : {
473 2221 : return _mproblem.numMatrixTags();
474 : }
475 :
476 : bool
477 5184 : DisplacedProblem::hasVariable(const std::string & var_name) const
478 : {
479 8446 : for (auto & nl : _displaced_solver_systems)
480 5184 : if (nl->hasVariable(var_name))
481 1922 : return true;
482 3262 : if (_displaced_aux->hasVariable(var_name))
483 3254 : return true;
484 :
485 8 : return false;
486 : }
487 :
488 : const MooseVariableFieldBase &
489 30624 : DisplacedProblem::getVariable(const THREAD_ID tid,
490 : const std::string & var_name,
491 : Moose::VarKindType expected_var_type,
492 : Moose::VarFieldType expected_var_field_type) const
493 : {
494 91872 : return getVariableHelper(tid,
495 : var_name,
496 : expected_var_type,
497 : expected_var_field_type,
498 30624 : _displaced_solver_systems,
499 61248 : *_displaced_aux);
500 : }
501 :
502 : MooseVariable &
503 42916 : DisplacedProblem::getStandardVariable(const THREAD_ID tid, const std::string & var_name)
504 : {
505 85678 : for (auto & nl : _displaced_solver_systems)
506 42916 : if (nl->hasVariable(var_name))
507 154 : return nl->getFieldVariable<Real>(tid, var_name);
508 42762 : if (_displaced_aux->hasVariable(var_name))
509 42762 : return _displaced_aux->getFieldVariable<Real>(tid, var_name);
510 :
511 0 : mooseError("No variable with name '" + var_name + "'");
512 : }
513 :
514 : MooseVariableFieldBase &
515 0 : DisplacedProblem::getActualFieldVariable(const THREAD_ID tid, const std::string & var_name)
516 : {
517 0 : for (auto & nl : _displaced_solver_systems)
518 0 : if (nl->hasVariable(var_name))
519 0 : return nl->getActualFieldVariable<Real>(tid, var_name);
520 0 : if (_displaced_aux->hasVariable(var_name))
521 0 : return _displaced_aux->getActualFieldVariable<Real>(tid, var_name);
522 :
523 0 : mooseError("No variable with name '" + var_name + "'");
524 : }
525 :
526 : VectorMooseVariable &
527 0 : DisplacedProblem::getVectorVariable(const THREAD_ID tid, const std::string & var_name)
528 : {
529 0 : for (auto & nl : _displaced_solver_systems)
530 0 : if (nl->hasVariable(var_name))
531 0 : return nl->getFieldVariable<RealVectorValue>(tid, var_name);
532 0 : if (_displaced_aux->hasVariable(var_name))
533 0 : return _displaced_aux->getFieldVariable<RealVectorValue>(tid, var_name);
534 :
535 0 : mooseError("No variable with name '" + var_name + "'");
536 : }
537 :
538 : ArrayMooseVariable &
539 0 : DisplacedProblem::getArrayVariable(const THREAD_ID tid, const std::string & var_name)
540 : {
541 0 : for (auto & nl : _displaced_solver_systems)
542 0 : if (nl->hasVariable(var_name))
543 0 : return nl->getFieldVariable<RealEigenVector>(tid, var_name);
544 0 : if (_displaced_aux->hasVariable(var_name))
545 0 : return _displaced_aux->getFieldVariable<RealEigenVector>(tid, var_name);
546 :
547 0 : mooseError("No variable with name '" + var_name + "'");
548 : }
549 :
550 : bool
551 2604 : DisplacedProblem::hasScalarVariable(const std::string & var_name) const
552 : {
553 5192 : for (auto & nl : _displaced_solver_systems)
554 2604 : if (nl->hasScalarVariable(var_name))
555 16 : return true;
556 2588 : if (_displaced_aux->hasScalarVariable(var_name))
557 0 : return true;
558 :
559 2588 : return false;
560 : }
561 :
562 : MooseVariableScalar &
563 16 : DisplacedProblem::getScalarVariable(const THREAD_ID tid, const std::string & var_name)
564 : {
565 16 : for (auto & nl : _displaced_solver_systems)
566 16 : if (nl->hasScalarVariable(var_name))
567 16 : return nl->getScalarVariable(tid, var_name);
568 0 : if (_displaced_aux->hasScalarVariable(var_name))
569 0 : return _displaced_aux->getScalarVariable(tid, var_name);
570 :
571 0 : mooseError("No variable with name '" + var_name + "'");
572 : }
573 :
574 : System &
575 12 : DisplacedProblem::getSystem(const std::string & var_name)
576 : {
577 24 : for (const auto sys_num : make_range(_eq.n_systems()))
578 : {
579 24 : auto & sys = _eq.get_system(sys_num);
580 24 : if (sys.has_variable(var_name))
581 12 : return sys;
582 : }
583 :
584 0 : mooseError("Unable to find a system containing the variable " + var_name);
585 : }
586 :
587 : void
588 3453 : DisplacedProblem::addVariable(const std::string & var_type,
589 : const std::string & name,
590 : InputParameters & parameters,
591 : const unsigned int nl_system_number)
592 : {
593 3453 : _displaced_solver_systems[nl_system_number]->addVariable(var_type, name, parameters);
594 3453 : }
595 :
596 : void
597 10166 : DisplacedProblem::addAuxVariable(const std::string & var_type,
598 : const std::string & name,
599 : InputParameters & parameters)
600 : {
601 10166 : _displaced_aux->addVariable(var_type, name, parameters);
602 10166 : }
603 :
604 : unsigned int
605 11463091 : DisplacedProblem::currentNlSysNum() const
606 : {
607 11463091 : return _mproblem.currentNlSysNum();
608 : }
609 :
610 : unsigned int
611 0 : DisplacedProblem::currentLinearSysNum() const
612 : {
613 0 : return _mproblem.currentLinearSysNum();
614 : }
615 :
616 : void
617 7563782 : DisplacedProblem::prepare(const Elem * elem, const THREAD_ID tid)
618 : {
619 15127540 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
620 : {
621 7563782 : _assembly[tid][nl_sys_num]->reinit(elem);
622 7563758 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
623 : // This method is called outside of residual/Jacobian callbacks during initial condition
624 : // evaluation
625 7563758 : if (!_mproblem.hasJacobian() || !_mproblem.constJacobian())
626 7563758 : _assembly[tid][nl_sys_num]->prepareJacobianBlock();
627 7563758 : _assembly[tid][nl_sys_num]->prepareResidual();
628 : }
629 :
630 7563758 : _displaced_aux->prepare(tid);
631 7563758 : }
632 :
633 : void
634 0 : DisplacedProblem::prepareNonlocal(const THREAD_ID tid)
635 : {
636 0 : _assembly[tid][currentNlSysNum()]->prepareNonlocal();
637 0 : }
638 :
639 : void
640 0 : DisplacedProblem::prepareFace(const Elem * /*elem*/, const THREAD_ID tid)
641 : {
642 0 : for (auto & nl : _displaced_solver_systems)
643 0 : nl->prepareFace(tid, true);
644 0 : _displaced_aux->prepareFace(tid, false);
645 0 : }
646 :
647 : void
648 0 : DisplacedProblem::prepare(const Elem * elem,
649 : unsigned int ivar,
650 : unsigned int jvar,
651 : const std::vector<dof_id_type> & dof_indices,
652 : const THREAD_ID tid)
653 : {
654 0 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
655 : {
656 0 : _assembly[tid][nl_sys_num]->reinit(elem);
657 0 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
658 : }
659 0 : _displaced_aux->prepare(tid);
660 0 : _assembly[tid][currentNlSysNum()]->prepareBlock(ivar, jvar, dof_indices);
661 0 : }
662 :
663 : void
664 31002 : DisplacedProblem::setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid)
665 : {
666 31002 : SubdomainID did = elem->subdomain_id();
667 62004 : for (auto & assembly : _assembly[tid])
668 31002 : assembly->setCurrentSubdomainID(did);
669 31002 : }
670 :
671 : void
672 129540 : DisplacedProblem::setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid)
673 : {
674 129540 : SubdomainID did = elem->neighbor_ptr(side)->subdomain_id();
675 259080 : for (auto & assembly : _assembly[tid])
676 129540 : assembly->setCurrentNeighborSubdomainID(did);
677 129540 : }
678 :
679 : void
680 0 : DisplacedProblem::prepareBlockNonlocal(unsigned int ivar,
681 : unsigned int jvar,
682 : const std::vector<dof_id_type> & idof_indices,
683 : const std::vector<dof_id_type> & jdof_indices,
684 : const THREAD_ID tid)
685 : {
686 0 : _assembly[tid][currentNlSysNum()]->prepareBlockNonlocal(ivar, jvar, idof_indices, jdof_indices);
687 0 : }
688 :
689 : void
690 51296 : DisplacedProblem::prepareAssembly(const THREAD_ID tid)
691 : {
692 51296 : _assembly[tid][currentNlSysNum()]->prepare();
693 51296 : }
694 :
695 : void
696 51296 : DisplacedProblem::prepareAssemblyNeighbor(const THREAD_ID tid)
697 : {
698 51296 : _assembly[tid][currentNlSysNum()]->prepareNeighbor();
699 51296 : }
700 :
701 : bool
702 4976 : DisplacedProblem::reinitDirac(const Elem * elem, const THREAD_ID tid)
703 : {
704 4976 : std::vector<Point> & points = _dirac_kernel_info.getPoints()[elem].first;
705 :
706 4976 : unsigned int n_points = points.size();
707 :
708 4976 : if (n_points)
709 : {
710 9952 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
711 : {
712 4976 : _assembly[tid][nl_sys_num]->reinitAtPhysical(elem, points);
713 4976 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
714 : }
715 4976 : _displaced_aux->prepare(tid);
716 :
717 4976 : reinitElem(elem, tid);
718 : }
719 :
720 4976 : _assembly[tid][currentNlSysNum()]->prepare();
721 :
722 4976 : return n_points > 0;
723 : }
724 :
725 : void
726 4431131 : DisplacedProblem::reinitElem(const Elem * elem, const THREAD_ID tid)
727 : {
728 8862262 : for (auto & nl : _displaced_solver_systems)
729 4431131 : nl->reinitElem(elem, tid);
730 4431131 : _displaced_aux->reinitElem(elem, tid);
731 4431131 : }
732 :
733 : void
734 0 : DisplacedProblem::reinitElemPhys(const Elem * elem,
735 : const std::vector<Point> & phys_points_in_elem,
736 : const THREAD_ID tid)
737 : {
738 : mooseAssert(_mesh.queryElemPtr(elem->id()) == elem,
739 : "Are you calling this method with a undisplaced mesh element?");
740 :
741 0 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
742 : {
743 0 : _assembly[tid][nl_sys_num]->reinitAtPhysical(elem, phys_points_in_elem);
744 0 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
745 0 : _assembly[tid][nl_sys_num]->prepare();
746 : }
747 0 : _displaced_aux->prepare(tid);
748 :
749 0 : reinitElem(elem, tid);
750 0 : }
751 :
752 : void
753 82248 : DisplacedProblem::reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid)
754 : {
755 164496 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
756 : {
757 82248 : _assembly[tid][nl_sys_num]->reinit(elem, side);
758 82248 : _displaced_solver_systems[nl_sys_num]->reinitElemFace(elem, side, tid);
759 : }
760 82248 : _displaced_aux->reinitElemFace(elem, side, tid);
761 82248 : }
762 :
763 : void
764 1001711 : DisplacedProblem::reinitNode(const Node * node, const THREAD_ID tid)
765 : {
766 2003422 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
767 : {
768 1001711 : _assembly[tid][nl_sys_num]->reinit(node);
769 1001711 : _displaced_solver_systems[nl_sys_num]->reinitNode(node, tid);
770 : }
771 1001711 : _displaced_aux->reinitNode(node, tid);
772 1001711 : }
773 :
774 : void
775 3388984 : DisplacedProblem::reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid)
776 : {
777 6777968 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
778 : {
779 3388984 : _assembly[tid][nl_sys_num]->reinit(node);
780 3388984 : _displaced_solver_systems[nl_sys_num]->reinitNodeFace(node, bnd_id, tid);
781 : }
782 3388984 : _displaced_aux->reinitNodeFace(node, bnd_id, tid);
783 3388984 : }
784 :
785 : void
786 0 : DisplacedProblem::reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid)
787 : {
788 0 : for (auto & nl : _displaced_solver_systems)
789 0 : nl->reinitNodes(nodes, tid);
790 0 : _displaced_aux->reinitNodes(nodes, tid);
791 0 : }
792 :
793 : void
794 0 : DisplacedProblem::reinitNodesNeighbor(const std::vector<dof_id_type> & nodes, const THREAD_ID tid)
795 : {
796 0 : for (auto & nl : _displaced_solver_systems)
797 0 : nl->reinitNodesNeighbor(nodes, tid);
798 0 : _displaced_aux->reinitNodesNeighbor(nodes, tid);
799 0 : }
800 :
801 : void
802 64740 : DisplacedProblem::reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid)
803 : {
804 64740 : reinitNeighbor(elem, side, tid, nullptr);
805 64740 : }
806 :
807 : void
808 129540 : DisplacedProblem::reinitNeighbor(const Elem * elem,
809 : unsigned int side,
810 : const THREAD_ID tid,
811 : const std::vector<Point> * neighbor_reference_points)
812 : {
813 129540 : setNeighborSubdomainID(elem, side, tid);
814 :
815 129540 : const Elem * neighbor = elem->neighbor_ptr(side);
816 129540 : unsigned int neighbor_side = neighbor->which_neighbor_am_i(elem);
817 :
818 259080 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
819 : {
820 129540 : _assembly[tid][nl_sys_num]->reinitElemAndNeighbor(
821 : elem, side, neighbor, neighbor_side, neighbor_reference_points);
822 129540 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
823 : // Called during stateful material property evaluation outside of solve
824 129540 : _assembly[tid][nl_sys_num]->prepareNeighbor();
825 : }
826 129540 : _displaced_aux->prepareNeighbor(tid);
827 :
828 259080 : for (auto & nl : _displaced_solver_systems)
829 : {
830 129540 : nl->reinitElemFace(elem, side, tid);
831 129540 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
832 : }
833 129540 : _displaced_aux->reinitElemFace(elem, side, tid);
834 129540 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
835 129540 : }
836 :
837 : void
838 51296 : DisplacedProblem::reinitNeighborPhys(const Elem * neighbor,
839 : unsigned int neighbor_side,
840 : const std::vector<Point> & physical_points,
841 : const THREAD_ID tid)
842 : {
843 : mooseAssert(_mesh.queryElemPtr(neighbor->id()) == neighbor,
844 : "Are you calling this method with a undisplaced mesh element?");
845 :
846 102592 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
847 : {
848 : // Reinit shape functions
849 51296 : _assembly[tid][nl_sys_num]->reinitNeighborAtPhysical(neighbor, neighbor_side, physical_points);
850 :
851 : // Set the neighbor dof indices
852 51296 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
853 : }
854 51296 : _displaced_aux->prepareNeighbor(tid);
855 :
856 51296 : prepareAssemblyNeighbor(tid);
857 :
858 : // Compute values at the points
859 102592 : for (auto & nl : _displaced_solver_systems)
860 51296 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
861 51296 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
862 51296 : }
863 :
864 : void
865 0 : DisplacedProblem::reinitNeighborPhys(const Elem * neighbor,
866 : const std::vector<Point> & physical_points,
867 : const THREAD_ID tid)
868 : {
869 : mooseAssert(_mesh.queryElemPtr(neighbor->id()) == neighbor,
870 : "Are you calling this method with a undisplaced mesh element?");
871 :
872 0 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
873 : {
874 : // Reinit shape functions
875 0 : _assembly[tid][nl_sys_num]->reinitNeighborAtPhysical(neighbor, physical_points);
876 :
877 : // Set the neighbor dof indices
878 0 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
879 : }
880 0 : _displaced_aux->prepareNeighbor(tid);
881 :
882 0 : prepareAssemblyNeighbor(tid);
883 :
884 : // Compute values at the points
885 0 : for (auto & nl : _displaced_solver_systems)
886 0 : nl->reinitNeighbor(neighbor, tid);
887 0 : _displaced_aux->reinitNeighbor(neighbor, tid);
888 0 : }
889 :
890 : void
891 64740 : DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem,
892 : unsigned int side,
893 : const THREAD_ID tid)
894 : {
895 64740 : reinitNeighbor(elem, side, tid);
896 :
897 64740 : const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
898 64740 : if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0)
899 480 : reinitLowerDElem(lower_d_elem, tid);
900 : else
901 : {
902 : // with mesh refinement, lower-dimensional element might be defined on neighbor side
903 64260 : auto & neighbor = _assembly[tid][currentNlSysNum()]->neighbor();
904 64260 : auto & neighbor_side = _assembly[tid][currentNlSysNum()]->neighborSide();
905 64260 : const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side);
906 64260 : if (lower_d_elem_neighbor &&
907 64260 : _mesh.interiorLowerDBlocks().count(lower_d_elem_neighbor->subdomain_id()) > 0)
908 : {
909 0 : auto qps = _assembly[tid][currentNlSysNum()]->qPointsFaceNeighbor().stdVector();
910 0 : std::vector<Point> reference_points;
911 0 : FEMap::inverse_map(
912 0 : lower_d_elem_neighbor->dim(), lower_d_elem_neighbor, qps, reference_points);
913 0 : reinitLowerDElem(lower_d_elem_neighbor, tid, &qps);
914 0 : }
915 : }
916 64740 : }
917 :
918 : void
919 113553 : DisplacedProblem::reinitScalars(const THREAD_ID tid,
920 : bool reinit_for_derivative_reordering /*=false*/)
921 : {
922 227106 : for (auto & nl : _displaced_solver_systems)
923 113553 : nl->reinitScalars(tid, reinit_for_derivative_reordering);
924 113553 : _displaced_aux->reinitScalars(tid, reinit_for_derivative_reordering);
925 113553 : }
926 :
927 : void
928 60 : DisplacedProblem::reinitOffDiagScalars(const THREAD_ID tid)
929 : {
930 60 : _assembly[tid][currentNlSysNum()]->prepareOffDiagScalar();
931 60 : }
932 :
933 : void
934 2359 : DisplacedProblem::getDiracElements(std::set<const Elem *> & elems)
935 : {
936 2359 : elems = _dirac_kernel_info.getElements();
937 2359 : }
938 :
939 : void
940 144140 : DisplacedProblem::clearDiracInfo()
941 : {
942 144140 : _dirac_kernel_info.clearPoints();
943 144140 : }
944 :
945 : void
946 4776 : DisplacedProblem::addResidual(const THREAD_ID tid)
947 : {
948 9552 : _assembly[tid][currentNlSysNum()]->addResidual(Assembly::GlobalDataKey{},
949 4776 : currentResidualVectorTags());
950 4776 : }
951 :
952 : void
953 61744 : DisplacedProblem::addResidualNeighbor(const THREAD_ID tid)
954 : {
955 123488 : _assembly[tid][currentNlSysNum()]->addResidualNeighbor(Assembly::GlobalDataKey{},
956 61744 : currentResidualVectorTags());
957 61744 : }
958 :
959 : void
960 61956 : DisplacedProblem::addResidualLower(const THREAD_ID tid)
961 : {
962 123912 : _assembly[tid][currentNlSysNum()]->addResidualLower(Assembly::GlobalDataKey{},
963 61956 : currentResidualVectorTags());
964 61956 : }
965 :
966 : void
967 35 : DisplacedProblem::addCachedResidualDirectly(NumericVector<Number> & residual, const THREAD_ID tid)
968 : {
969 70 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
970 35 : _displaced_solver_systems[currentNlSysNum()]->timeVectorTag()))
971 0 : _assembly[tid][currentNlSysNum()]->addCachedResidualDirectly(
972 : residual,
973 0 : Assembly::GlobalDataKey{},
974 0 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->timeVectorTag()));
975 :
976 70 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
977 35 : _displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag()))
978 105 : _assembly[tid][currentNlSysNum()]->addCachedResidualDirectly(
979 : residual,
980 70 : Assembly::GlobalDataKey{},
981 35 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag()));
982 :
983 35 : std::vector<VectorTag> extra_residual_vector_tags;
984 35 : extra_residual_vector_tags.reserve(currentResidualVectorTags().size());
985 35 : const auto time_tag = _displaced_solver_systems[currentNlSysNum()]->timeVectorTag();
986 35 : const auto non_time_tag = _displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag();
987 105 : for (const auto & vector_tag : currentResidualVectorTags())
988 70 : if (vector_tag._id != time_tag && vector_tag._id != non_time_tag)
989 35 : extra_residual_vector_tags.push_back(vector_tag);
990 :
991 : // Flush extra vector tag caches (e.g. from extra_vector_tags on NodalConstraints)
992 : // to their respective system vectors after the standard TIME/NONTIME caches above.
993 : // Without this, NodalConstraint contributions to extra vector tags are silently
994 : // discarded by the blanket clearCachedResiduals.
995 35 : _assembly[tid][currentNlSysNum()]->addCachedResiduals(Assembly::GlobalDataKey{},
996 : extra_residual_vector_tags);
997 :
998 : // We do this because by adding the cached residual directly, we cannot ensure that all of the
999 : // cached residuals are emptied after only the two add calls above
1000 35 : _assembly[tid][currentNlSysNum()]->clearCachedResiduals(Assembly::GlobalDataKey{});
1001 35 : }
1002 :
1003 : void
1004 0 : DisplacedProblem::setResidual(NumericVector<Number> & residual, const THREAD_ID tid)
1005 : {
1006 0 : _assembly[tid][currentNlSysNum()]->setResidual(
1007 : residual,
1008 0 : Assembly::GlobalDataKey{},
1009 0 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->residualVectorTag()));
1010 0 : }
1011 :
1012 : void
1013 0 : DisplacedProblem::setResidualNeighbor(NumericVector<Number> & residual, const THREAD_ID tid)
1014 : {
1015 0 : _assembly[tid][currentNlSysNum()]->setResidualNeighbor(
1016 : residual,
1017 0 : Assembly::GlobalDataKey{},
1018 0 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->residualVectorTag()));
1019 0 : }
1020 :
1021 : void
1022 200 : DisplacedProblem::addJacobian(const THREAD_ID tid)
1023 : {
1024 200 : _assembly[tid][currentNlSysNum()]->addJacobian(Assembly::GlobalDataKey{});
1025 200 : }
1026 :
1027 : void
1028 0 : DisplacedProblem::addJacobianNonlocal(const THREAD_ID tid)
1029 : {
1030 0 : _assembly[tid][currentNlSysNum()]->addJacobianNonlocal(Assembly::GlobalDataKey{});
1031 0 : }
1032 :
1033 : void
1034 44 : DisplacedProblem::addJacobianNeighbor(const THREAD_ID tid)
1035 : {
1036 44 : _assembly[tid][currentNlSysNum()]->addJacobianNeighbor(Assembly::GlobalDataKey{});
1037 44 : }
1038 :
1039 : void
1040 3072 : DisplacedProblem::addJacobianNeighborLowerD(const THREAD_ID tid)
1041 : {
1042 3072 : _assembly[tid][currentNlSysNum()]->addJacobianNeighborLowerD(Assembly::GlobalDataKey{});
1043 3072 : }
1044 :
1045 : void
1046 192 : DisplacedProblem::addJacobianLowerD(const THREAD_ID tid)
1047 : {
1048 192 : _assembly[tid][currentNlSysNum()]->addJacobianLowerD(Assembly::GlobalDataKey{});
1049 192 : }
1050 :
1051 : void
1052 0 : DisplacedProblem::cacheJacobianNonlocal(const THREAD_ID tid)
1053 : {
1054 0 : _assembly[tid][currentNlSysNum()]->cacheJacobianNonlocal(Assembly::GlobalDataKey{});
1055 0 : }
1056 :
1057 : void
1058 0 : DisplacedProblem::addJacobianBlockTags(SparseMatrix<Number> & jacobian,
1059 : unsigned int ivar,
1060 : unsigned int jvar,
1061 : const DofMap & dof_map,
1062 : std::vector<dof_id_type> & dof_indices,
1063 : const std::set<TagID> & tags,
1064 : const THREAD_ID tid)
1065 : {
1066 0 : _assembly[tid][currentNlSysNum()]->addJacobianBlockTags(
1067 0 : jacobian, ivar, jvar, dof_map, dof_indices, Assembly::GlobalDataKey{}, tags);
1068 0 : }
1069 :
1070 : void
1071 0 : DisplacedProblem::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
1072 : unsigned int ivar,
1073 : unsigned int jvar,
1074 : const DofMap & dof_map,
1075 : const std::vector<dof_id_type> & idof_indices,
1076 : const std::vector<dof_id_type> & jdof_indices,
1077 : const std::set<TagID> & tags,
1078 : const THREAD_ID tid)
1079 : {
1080 0 : _assembly[tid][currentNlSysNum()]->addJacobianBlockNonlocalTags(
1081 0 : jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices, Assembly::GlobalDataKey{}, tags);
1082 0 : }
1083 :
1084 : void
1085 0 : DisplacedProblem::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
1086 : unsigned int ivar,
1087 : unsigned int jvar,
1088 : const DofMap & dof_map,
1089 : std::vector<dof_id_type> & dof_indices,
1090 : std::vector<dof_id_type> & neighbor_dof_indices,
1091 : const std::set<TagID> & tags,
1092 : const THREAD_ID tid)
1093 : {
1094 0 : _assembly[tid][currentNlSysNum()]->addJacobianNeighborTags(jacobian,
1095 : ivar,
1096 : jvar,
1097 : dof_map,
1098 : dof_indices,
1099 : neighbor_dof_indices,
1100 0 : Assembly::GlobalDataKey{},
1101 : tags);
1102 0 : }
1103 :
1104 : void
1105 745429 : DisplacedProblem::prepareShapes(unsigned int var, const THREAD_ID tid)
1106 : {
1107 745429 : _assembly[tid][currentNlSysNum()]->copyShapes(var);
1108 745429 : }
1109 :
1110 : void
1111 8758 : DisplacedProblem::prepareFaceShapes(unsigned int var, const THREAD_ID tid)
1112 : {
1113 8758 : _assembly[tid][currentNlSysNum()]->copyFaceShapes(var);
1114 8758 : }
1115 :
1116 : void
1117 3352 : DisplacedProblem::prepareNeighborShapes(unsigned int var, const THREAD_ID tid)
1118 : {
1119 3352 : _assembly[tid][currentNlSysNum()]->copyNeighborShapes(var);
1120 3352 : }
1121 :
1122 : void
1123 4209 : DisplacedProblem::updateGeomSearch(GeometricSearchData::GeometricSearchType type)
1124 : {
1125 21045 : TIME_SECTION("updateGeometricSearch", 3, "Updating Displaced GeometricSearch");
1126 :
1127 4209 : _geometric_search_data.update(type);
1128 4209 : }
1129 :
1130 : void
1131 569 : DisplacedProblem::meshChanged(const bool contract_mesh, const bool clean_refinement_flags)
1132 : {
1133 : // The mesh changed. The displaced equations system object only holds Systems, so calling
1134 : // EquationSystems::reinit only prolongs/restricts the solution vectors, which is something that
1135 : // needs to happen for every step of mesh adaptivity.
1136 569 : _eq.reinit();
1137 569 : if (contract_mesh)
1138 : // Once vectors are restricted, we can delete children of coarsened elements
1139 444 : _mesh.getMesh().contract();
1140 569 : if (clean_refinement_flags)
1141 : {
1142 : // Finally clean refinement flags so that if someone tries to project vectors again without
1143 : // an intervening mesh refinement to clean flags they won't run into trouble
1144 444 : MeshRefinement refinement(_mesh.getMesh());
1145 444 : refinement.clean_refinement_flags();
1146 444 : }
1147 :
1148 : // Since the mesh has changed, we need to make sure that we update any of our
1149 : // MOOSE-system specific data.
1150 1138 : for (auto & nl : _displaced_solver_systems)
1151 569 : nl->reinit();
1152 569 : _displaced_aux->reinit();
1153 :
1154 : // We've performed some mesh adaptivity. We need to
1155 : // clear any quadrature nodes such that when we build the boundary node lists in
1156 : // MooseMesh::meshChanged we don't have any extraneous extra boundary nodes lying around
1157 569 : _mesh.clearQuadratureNodes();
1158 :
1159 569 : _mesh.meshChanged();
1160 :
1161 : // Before performing mesh adaptivity we un-displaced the mesh. We need to re-displace the mesh and
1162 : // then reinitialize GeometricSearchData such that we have all the correct geometric information
1163 : // for the changed mesh
1164 569 : updateMesh(/*mesh_changing=*/true);
1165 569 : }
1166 :
1167 : void
1168 808450 : DisplacedProblem::addGhostedElem(dof_id_type elem_id)
1169 : {
1170 808450 : _mproblem.addGhostedElem(elem_id);
1171 808450 : }
1172 :
1173 : void
1174 26702 : DisplacedProblem::addGhostedBoundary(BoundaryID boundary_id)
1175 : {
1176 26702 : _mproblem.addGhostedBoundary(boundary_id);
1177 26702 : }
1178 :
1179 : void
1180 0 : DisplacedProblem::ghostGhostedBoundaries()
1181 : {
1182 0 : _mproblem.ghostGhostedBoundaries();
1183 0 : }
1184 :
1185 : MooseMesh &
1186 408883 : DisplacedProblem::refMesh()
1187 : {
1188 408883 : return _ref_mesh;
1189 : }
1190 :
1191 : bool
1192 0 : DisplacedProblem::solverSystemConverged(const unsigned int sys_num)
1193 : {
1194 0 : return _mproblem.converged(sys_num);
1195 : }
1196 :
1197 : bool
1198 0 : DisplacedProblem::computingPreSMOResidual(const unsigned int nl_sys_num) const
1199 : {
1200 0 : return _mproblem.computingPreSMOResidual(nl_sys_num);
1201 : }
1202 :
1203 : void
1204 0 : DisplacedProblem::onTimestepBegin()
1205 : {
1206 0 : }
1207 :
1208 : void
1209 0 : DisplacedProblem::onTimestepEnd()
1210 : {
1211 0 : }
1212 :
1213 : void
1214 567 : DisplacedProblem::undisplaceMesh()
1215 : {
1216 : // If undisplaceMesh() is called during initial adaptivity, it is
1217 : // not valid to call _mesh.getActiveSemiLocalNodeRange() since it is
1218 : // not set up yet. So we are creating the Range by hand.
1219 : //
1220 : // We must undisplace *all* our nodes to the _ref_mesh
1221 : // configuration, not just the local ones, since the partitioners
1222 : // require this. We are using the GRAIN_SIZE=1 from MooseMesh.C,
1223 : // not sure how this value was decided upon.
1224 : //
1225 : // (DRG: The grainsize parameter is ultimately passed to TBB to help
1226 : // it choose how to split up the range. A grainsize of 1 says "split
1227 : // it as much as you want". Years ago I experimentally found that it
1228 : // didn't matter much and that using 1 was fine.)
1229 : //
1230 : // Note: we don't have to invalidate/update as much stuff as
1231 : // DisplacedProblem::updateMesh() does, since this will be handled
1232 : // by a later call to updateMesh().
1233 1134 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
1234 1134 : _mesh.getMesh().nodes_end(),
1235 567 : /*grainsize=*/1);
1236 :
1237 567 : ResetDisplacedMeshThread rdmt(_mproblem, *this);
1238 :
1239 : // Undisplace the mesh using threads.
1240 567 : Threads::parallel_reduce(node_range, rdmt);
1241 567 : }
1242 :
1243 : LineSearch *
1244 0 : DisplacedProblem::getLineSearch()
1245 : {
1246 0 : return _mproblem.getLineSearch();
1247 : }
1248 :
1249 : const CouplingMatrix *
1250 0 : DisplacedProblem::couplingMatrix(const unsigned int nl_sys_num) const
1251 : {
1252 0 : return _mproblem.couplingMatrix(nl_sys_num);
1253 : }
1254 :
1255 : bool
1256 502089 : DisplacedProblem::computingScalingJacobian() const
1257 : {
1258 502089 : return _mproblem.computingScalingJacobian();
1259 : }
1260 :
1261 : bool
1262 0 : DisplacedProblem::computingScalingResidual() const
1263 : {
1264 0 : return _mproblem.computingScalingResidual();
1265 : }
1266 :
1267 : void
1268 2013 : DisplacedProblem::initialSetup()
1269 : {
1270 2013 : SubProblem::initialSetup();
1271 :
1272 4026 : for (auto & nl : _displaced_solver_systems)
1273 2013 : nl->initialSetup();
1274 2013 : _displaced_aux->initialSetup();
1275 2013 : }
1276 :
1277 : void
1278 30674 : DisplacedProblem::timestepSetup()
1279 : {
1280 30674 : SubProblem::timestepSetup();
1281 :
1282 61348 : for (auto & nl : _displaced_solver_systems)
1283 30674 : nl->timestepSetup();
1284 30674 : _displaced_aux->timestepSetup();
1285 30674 : }
1286 :
1287 : void
1288 145157 : DisplacedProblem::customSetup(const ExecFlagType & exec_type)
1289 : {
1290 145157 : SubProblem::customSetup(exec_type);
1291 :
1292 290314 : for (auto & nl : _displaced_solver_systems)
1293 145157 : nl->customSetup(exec_type);
1294 145157 : _displaced_aux->customSetup(exec_type);
1295 145157 : }
1296 :
1297 : void
1298 124149 : DisplacedProblem::residualSetup()
1299 : {
1300 124149 : SubProblem::residualSetup();
1301 :
1302 248298 : for (auto & nl : _displaced_solver_systems)
1303 124149 : nl->residualSetup();
1304 124149 : _displaced_aux->residualSetup();
1305 124149 : }
1306 :
1307 : void
1308 21144 : DisplacedProblem::jacobianSetup()
1309 : {
1310 21144 : SubProblem::jacobianSetup();
1311 :
1312 42288 : for (auto & nl : _displaced_solver_systems)
1313 21144 : nl->jacobianSetup();
1314 21144 : _displaced_aux->jacobianSetup();
1315 21144 : }
1316 :
1317 : void
1318 284 : DisplacedProblem::haveADObjects(const bool have_ad_objects)
1319 : {
1320 284 : _have_ad_objects = have_ad_objects;
1321 284 : _mproblem.SubProblem::haveADObjects(have_ad_objects);
1322 284 : }
1323 :
1324 : std::pair<bool, unsigned int>
1325 30624 : DisplacedProblem::determineSolverSystem(const std::string & var_name,
1326 : const bool error_if_not_found) const
1327 : {
1328 30624 : return _mproblem.determineSolverSystem(var_name, error_if_not_found);
1329 : }
1330 :
1331 : Assembly &
1332 52741260 : DisplacedProblem::assembly(const THREAD_ID tid, const unsigned int sys_num)
1333 : {
1334 : mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
1335 : mooseAssert(sys_num < _assembly[tid].size(),
1336 : "System number larger than the assembly container size");
1337 52741260 : return *_assembly[tid][sys_num];
1338 : }
1339 :
1340 : const Assembly &
1341 44805 : DisplacedProblem::assembly(const THREAD_ID tid, const unsigned int sys_num) const
1342 : {
1343 : mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
1344 : mooseAssert(sys_num < _assembly[tid].size(),
1345 : "System number larger than the assembly container size");
1346 44805 : return *_assembly[tid][sys_num];
1347 : }
1348 :
1349 : std::size_t
1350 9179383 : DisplacedProblem::numNonlinearSystems() const
1351 : {
1352 9179383 : return _mproblem.numNonlinearSystems();
1353 : }
1354 :
1355 : std::size_t
1356 0 : DisplacedProblem::numLinearSystems() const
1357 : {
1358 0 : return _mproblem.numLinearSystems();
1359 : }
1360 :
1361 : std::size_t
1362 359320 : DisplacedProblem::numSolverSystems() const
1363 : {
1364 359320 : return _mproblem.numSolverSystems();
1365 : }
1366 :
1367 : const std::vector<VectorTag> &
1368 8214459 : DisplacedProblem::currentResidualVectorTags() const
1369 : {
1370 8214459 : return _mproblem.currentResidualVectorTags();
1371 : }
1372 :
1373 : bool
1374 45474904 : DisplacedProblem::safeAccessTaggedMatrices() const
1375 : {
1376 45474904 : return _mproblem.safeAccessTaggedMatrices();
1377 : }
1378 :
1379 : bool
1380 343 : DisplacedProblem::safeAccessTaggedVectors() const
1381 : {
1382 343 : return _mproblem.safeAccessTaggedVectors();
1383 : }
1384 :
1385 : void
1386 0 : DisplacedProblem::needFV()
1387 : {
1388 0 : _mproblem.needFV();
1389 0 : }
1390 :
1391 : bool
1392 295531 : DisplacedProblem::haveFV() const
1393 : {
1394 295531 : return _mproblem.haveFV();
1395 : }
1396 :
1397 : bool
1398 1544138 : DisplacedProblem::hasNonlocalCoupling() const
1399 : {
1400 1544138 : return _mproblem.hasNonlocalCoupling();
1401 : }
1402 :
1403 : unsigned int
1404 0 : DisplacedProblem::nlSysNum(const NonlinearSystemName & nl_sys_name) const
1405 : {
1406 0 : return _mproblem.nlSysNum(nl_sys_name);
1407 : }
1408 :
1409 : unsigned int
1410 0 : DisplacedProblem::linearSysNum(const LinearSystemName & sys_name) const
1411 : {
1412 0 : return _mproblem.linearSysNum(sys_name);
1413 : }
1414 :
1415 : unsigned int
1416 0 : DisplacedProblem::solverSysNum(const SolverSystemName & sys_name) const
1417 : {
1418 0 : return _mproblem.solverSysNum(sys_name);
1419 : }
1420 :
1421 : const libMesh::CouplingMatrix &
1422 2222 : DisplacedProblem::nonlocalCouplingMatrix(const unsigned i) const
1423 : {
1424 2222 : return _mproblem.nonlocalCouplingMatrix(i);
1425 : }
1426 :
1427 : bool
1428 0 : DisplacedProblem::checkNonlocalCouplingRequirement() const
1429 : {
1430 0 : return _mproblem.checkNonlocalCouplingRequirement();
1431 : }
1432 :
1433 148216 : DisplacedProblem::UpdateDisplacedMeshThread::UpdateDisplacedMeshThread(
1434 148216 : FEProblemBase & fe_problem, DisplacedProblem & displaced_problem)
1435 : : ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>(fe_problem),
1436 148216 : _displaced_problem(displaced_problem),
1437 296432 : _ref_mesh(_displaced_problem.refMesh()),
1438 148216 : _nl_soln(_displaced_problem._nl_solution),
1439 148216 : _aux_soln(*_displaced_problem._aux_solution),
1440 148216 : _has_displacement(false)
1441 : {
1442 148216 : this->init();
1443 148216 : }
1444 :
1445 14503 : DisplacedProblem::UpdateDisplacedMeshThread::UpdateDisplacedMeshThread(
1446 14503 : UpdateDisplacedMeshThread & x, Threads::split split)
1447 : : ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>(x, split),
1448 14503 : _displaced_problem(x._displaced_problem),
1449 14503 : _ref_mesh(x._ref_mesh),
1450 14503 : _nl_soln(x._nl_soln),
1451 14503 : _aux_soln(x._aux_soln),
1452 14503 : _sys_to_nonghost_and_ghost_soln(x._sys_to_nonghost_and_ghost_soln),
1453 14503 : _sys_to_var_num_and_direction(x._sys_to_var_num_and_direction),
1454 14503 : _has_displacement(x._has_displacement)
1455 : {
1456 14503 : }
1457 :
1458 : void
1459 148216 : DisplacedProblem::UpdateDisplacedMeshThread::init()
1460 : {
1461 148216 : std::vector<std::string> & displacement_variables = _displaced_problem._displacements;
1462 148216 : unsigned int num_displacements = displacement_variables.size();
1463 148216 : auto & es = _displaced_problem.es();
1464 :
1465 148216 : _sys_to_var_num_and_direction.clear();
1466 148216 : _sys_to_nonghost_and_ghost_soln.clear();
1467 :
1468 495580 : for (unsigned int i = 0; i < num_displacements; i++)
1469 : {
1470 347364 : std::string displacement_name = displacement_variables[i];
1471 :
1472 456078 : for (const auto sys_num : make_range(es.n_systems()))
1473 : {
1474 456078 : auto & sys = es.get_system(sys_num);
1475 456078 : if (sys.has_variable(displacement_name))
1476 : {
1477 347364 : auto & val = _sys_to_var_num_and_direction[sys.number()];
1478 347364 : val.first.push_back(sys.variable_number(displacement_name));
1479 347364 : val.second.push_back(i);
1480 347364 : break;
1481 : }
1482 : }
1483 347364 : }
1484 :
1485 298550 : for (const auto & pr : _sys_to_var_num_and_direction)
1486 : {
1487 150334 : auto & sys = es.get_system(pr.first);
1488 : mooseAssert(sys.number() <= _nl_soln.size(),
1489 : "The system number should always be less than or equal to the number of nonlinear "
1490 : "systems. If it is equal, then this system is the auxiliary system");
1491 : const NumericVector<Number> * const nonghost_soln =
1492 150334 : sys.number() < _nl_soln.size() ? _nl_soln[sys.number()] : &_aux_soln;
1493 300668 : _sys_to_nonghost_and_ghost_soln.emplace(
1494 150334 : sys.number(),
1495 0 : std::make_pair(nonghost_soln,
1496 300668 : NumericVector<Number>::build(nonghost_soln->comm()).release()));
1497 : }
1498 :
1499 148216 : ConstNodeRange node_range(_ref_mesh.getMesh().nodes_begin(), _ref_mesh.getMesh().nodes_end());
1500 :
1501 298550 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1502 : {
1503 150334 : auto & sys = es.get_system(sys_num);
1504 : AllNodesSendListThread send_list(
1505 150334 : this->_fe_problem, _ref_mesh, var_num_and_direction.first, sys);
1506 150334 : Threads::parallel_reduce(node_range, send_list);
1507 150334 : send_list.unique();
1508 150334 : auto & [soln, ghost_soln] = libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num);
1509 300668 : ghost_soln->init(
1510 150334 : soln->size(), soln->local_size(), send_list.send_list(), true, libMesh::GHOSTED);
1511 150334 : soln->localize(*ghost_soln, send_list.send_list());
1512 150334 : }
1513 :
1514 148216 : _has_displacement = false;
1515 148216 : }
1516 :
1517 : void
1518 18026081 : DisplacedProblem::UpdateDisplacedMeshThread::onNode(NodeRange::const_iterator & nd)
1519 : {
1520 18026081 : Node & displaced_node = *(*nd);
1521 :
1522 18026081 : Node & reference_node = _ref_mesh.nodeRef(displaced_node.id());
1523 :
1524 36472836 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1525 : {
1526 18446755 : auto & var_numbers = var_num_and_direction.first;
1527 18446755 : auto & directions = var_num_and_direction.second;
1528 56813405 : for (const auto i : index_range(var_numbers))
1529 : {
1530 38366650 : const auto direction = directions[i];
1531 38366650 : if (reference_node.n_dofs(sys_num, var_numbers[i]) > 0)
1532 : {
1533 38366522 : Real coord = reference_node(direction) +
1534 76733044 : (*libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num).second)(
1535 38366522 : reference_node.dof_number(sys_num, var_numbers[i], 0));
1536 38366522 : if (displaced_node(direction) != coord)
1537 : {
1538 3124798 : displaced_node(direction) = coord;
1539 3124798 : _has_displacement = true;
1540 : }
1541 : }
1542 : }
1543 : }
1544 18026081 : }
|
