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 18731 : DisplacedProblem::validParams()
35 : {
36 18731 : InputParameters params = SubProblem::validParams();
37 18731 : params.addClassDescription(
38 : "A Problem object for providing access to the displaced finite element "
39 : "mesh and associated variables.");
40 18731 : params.addPrivateParam<MooseMesh *>("mesh");
41 18731 : params.addPrivateParam<std::vector<std::string>>("displacements", {});
42 18731 : return params;
43 0 : }
44 :
45 2233 : DisplacedProblem::DisplacedProblem(const InputParameters & parameters)
46 : : SubProblem(parameters),
47 2233 : _mproblem(parameters.have_parameter<FEProblemBase *>("_fe_problem_base")
48 4466 : ? *getParam<FEProblemBase *>("_fe_problem_base")
49 2233 : : *getParam<FEProblem *>("_fe_problem")),
50 2233 : _mesh(*getParam<MooseMesh *>("mesh")),
51 2233 : _eq(_mesh),
52 2233 : _ref_mesh(_mproblem.mesh()),
53 2233 : _displacements(getParam<std::vector<std::string>>("displacements")),
54 4466 : _geometric_search_data(*this, _mesh)
55 :
56 : {
57 : // Disable refinement/coarsening in EquationSystems::reinit because we already do this ourselves
58 2233 : _eq.disable_refine_in_reinit();
59 :
60 : // TODO: Move newAssemblyArray further up to SubProblem so that we can use it here
61 2233 : unsigned int n_threads = libMesh::n_threads();
62 :
63 2233 : _assembly.resize(n_threads);
64 4466 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
65 : {
66 2233 : _displaced_solver_systems.emplace_back(std::make_unique<DisplacedSystem>(
67 : *this,
68 : _mproblem,
69 2233 : _mproblem.getNonlinearSystemBase(nl_sys_num),
70 4466 : "displaced_" + _mproblem.getNonlinearSystemBase(nl_sys_num).name() + "_" +
71 2233 : std::to_string(nl_sys_num),
72 2233 : Moose::VAR_SOLVER));
73 2233 : auto & displaced_nl = _displaced_solver_systems.back();
74 :
75 4665 : for (unsigned int i = 0; i < n_threads; ++i)
76 2432 : _assembly[i].emplace_back(std::make_unique<Assembly>(*displaced_nl, i));
77 : }
78 :
79 2233 : _nl_solution.resize(_displaced_solver_systems.size(), nullptr);
80 :
81 : _displaced_aux =
82 4466 : std::make_unique<DisplacedSystem>(*this,
83 : _mproblem,
84 2233 : _mproblem.getAuxiliarySystem(),
85 4466 : "displaced_" + _mproblem.getAuxiliarySystem().name(),
86 4466 : 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 2233 : automaticScaling(_mproblem.automaticScaling());
100 :
101 2233 : _mesh.setCoordData(_ref_mesh);
102 2233 : }
103 :
104 4386 : DisplacedProblem::~DisplacedProblem() = default;
105 :
106 : bool
107 124820527 : DisplacedProblem::isTransient() const
108 : {
109 124820527 : return _mproblem.isTransient();
110 : }
111 :
112 : std::set<dof_id_type> &
113 9442 : DisplacedProblem::ghostedElems()
114 : {
115 9442 : return _mproblem.ghostedElems();
116 : }
117 :
118 : void
119 2233 : 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 4665 : for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
127 4864 : for (const auto sys_num : index_range(_assembly[tid]))
128 2432 : _assembly[tid][sys_num]->createQRules(
129 : type, order, volume_order, face_order, block, allow_negative_qweights);
130 2233 : }
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 2233 : DisplacedProblem::init()
150 : {
151 4665 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
152 : {
153 4864 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
154 2432 : _assembly[tid][nl_sys_num]->init(_mproblem.couplingMatrix(nl_sys_num));
155 :
156 4864 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
157 : {
158 2432 : std::vector<std::pair<unsigned int, unsigned short>> disp_numbers_and_directions;
159 7680 : for (const auto direction : index_range(_displacements))
160 : {
161 5248 : const auto & disp_string = _displacements[direction];
162 5248 : const auto & disp_variable = getVariable(tid, disp_string);
163 5248 : if (disp_variable.sys().number() == nl_sys_num)
164 2352 : disp_numbers_and_directions.push_back(
165 4704 : std::make_pair(disp_variable.number(), cast_int<unsigned short>(direction)));
166 : }
167 2432 : _assembly[tid][nl_sys_num]->assignDisplacements(std::move(disp_numbers_and_directions));
168 2432 : }
169 : }
170 :
171 4466 : for (auto & nl : _displaced_solver_systems)
172 : {
173 2233 : nl->dofMap().attach_extra_send_list_function(&extraSendList, nl.get());
174 2233 : nl->preInit();
175 : }
176 :
177 2233 : _displaced_aux->dofMap().attach_extra_send_list_function(&extraSendList, _displaced_aux.get());
178 2233 : _displaced_aux->preInit();
179 :
180 : {
181 2233 : TIME_SECTION("eq::init", 2, "Initializing Displaced Equation System");
182 2233 : _eq.init();
183 2233 : }
184 :
185 4466 : for (auto & nl : _displaced_solver_systems)
186 2233 : nl->postInit();
187 2233 : _displaced_aux->postInit();
188 :
189 2233 : _mesh.meshChanged();
190 :
191 2233 : if (haveFV())
192 13 : _mesh.setupFiniteVolumeMeshData();
193 2233 : }
194 :
195 : void
196 142 : DisplacedProblem::initAdaptivity()
197 : {
198 142 : }
199 :
200 : void
201 1895 : DisplacedProblem::addTimeIntegrator()
202 : {
203 3790 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
204 1895 : _displaced_solver_systems[nl_sys_num]->copyTimeIntegrators(
205 1895 : _mproblem.getNonlinearSystemBase(nl_sys_num));
206 1895 : _displaced_aux->copyTimeIntegrators(_mproblem.getAuxiliarySystem());
207 1895 : }
208 :
209 : void
210 11 : DisplacedProblem::saveOldSolutions()
211 : {
212 22 : for (auto & displaced_nl : _displaced_solver_systems)
213 11 : displaced_nl->saveOldSolutions();
214 11 : _displaced_aux->saveOldSolutions();
215 11 : }
216 :
217 : void
218 11 : DisplacedProblem::restoreOldSolutions()
219 : {
220 22 : for (auto & displaced_nl : _displaced_solver_systems)
221 11 : displaced_nl->restoreOldSolutions();
222 11 : _displaced_aux->restoreOldSolutions();
223 11 : }
224 :
225 : void
226 535340 : DisplacedProblem::syncAuxSolution(const NumericVector<Number> & aux_soln)
227 : {
228 535340 : (*_displaced_aux->sys().solution) = aux_soln;
229 535340 : _displaced_aux->update();
230 535340 : }
231 :
232 : void
233 338980 : DisplacedProblem::syncSolutions()
234 : {
235 338980 : TIME_SECTION("syncSolutions", 5, "Syncing Displaced Solutions");
236 :
237 677960 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
238 : {
239 338980 : 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 338980 : (*displaced_nl->sys().solution) =
244 677960 : *_mproblem.getNonlinearSystemBase(displaced_nl->number()).currentSolution();
245 338980 : displaced_nl->update();
246 : }
247 338980 : syncAuxSolution(*_mproblem.getAuxiliarySystem().currentSolution());
248 338980 : }
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 163565 : DisplacedProblem::updateMesh(bool mesh_changing)
267 : {
268 163565 : 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 163565 : if (!mesh_changing)
275 163008 : syncSolutions();
276 :
277 327130 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
278 163565 : _nl_solution[nl_sys_num] = _displaced_solver_systems[nl_sys_num]->sys().solution.get();
279 163565 : _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 163565 : if (_mesh.getMesh().is_serial() && !this->refMesh().getMesh().is_serial())
287 0 : this->refMesh().getMesh().allgather();
288 :
289 163565 : if (_mesh.getMesh().is_serial_on_zero() && !this->refMesh().getMesh().is_serial_on_zero())
290 0 : this->refMesh().getMesh().gather_to_zero();
291 :
292 163565 : UpdateDisplacedMeshThread udmt(_mproblem, *this);
293 :
294 : // We displace all nodes, not just semilocal nodes, because
295 : // parallel-inconsistent mesh geometry makes libMesh cry.
296 327130 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
297 327130 : _mesh.getMesh().nodes_end(),
298 163565 : /*grainsize=*/1);
299 :
300 163565 : Threads::parallel_reduce(node_range, udmt);
301 : // Displacement of the mesh has invalidated the point locator data (e.g. bounding boxes)
302 163565 : _mesh.getMesh().clear_point_locator();
303 :
304 : // The mesh has changed. Face information normals, areas, etc. must be re-calculated
305 163565 : if (haveFV())
306 13 : _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 163565 : if (mesh_changing)
317 557 : _geometric_search_data.reinit();
318 : else
319 163008 : _geometric_search_data.update();
320 : }
321 0 : catch (MooseException & e)
322 : {
323 0 : _mproblem.setException(e.what());
324 0 : }
325 :
326 163561 : if (udmt.hasDisplacement())
327 70072 : _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 163561 : _mproblem.checkExceptionAndStopSolve(/*print_message=*/false);
333 :
334 : // Since the Mesh changed, update the PointLocator object used by DiracKernels.
335 163561 : _dirac_kernel_info.updatePointLocator(_mesh);
336 163561 : }
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 6187173 : DisplacedProblem::getVectorTag(const TagID tag_id) const
394 : {
395 6187173 : return _mproblem.getVectorTag(tag_id);
396 : }
397 :
398 : TagID
399 115335 : DisplacedProblem::getVectorTagID(const TagName & tag_name) const
400 : {
401 115335 : 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 237 : DisplacedProblem::vectorTagExists(const TagID tag_id) const
412 : {
413 237 : 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 48391 : DisplacedProblem::numVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
424 : {
425 48391 : return _mproblem.numVectorTags(type);
426 : }
427 :
428 : const std::vector<VectorTag> &
429 2432 : DisplacedProblem::getVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
430 : {
431 2432 : return _mproblem.getVectorTags(type);
432 : }
433 :
434 : Moose::VectorTagType
435 128305016 : DisplacedProblem::vectorTagType(const TagID tag_id) const
436 : {
437 128305016 : 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 1613 : DisplacedProblem::getMatrixTagID(const TagName & tag_name) const
448 : {
449 1613 : 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 388 : DisplacedProblem::matrixTagExists(TagID tag_id) const
466 : {
467 388 : return _mproblem.matrixTagExists(tag_id);
468 : }
469 :
470 : unsigned int
471 50823 : DisplacedProblem::numMatrixTags() const
472 : {
473 50823 : return _mproblem.numMatrixTags();
474 : }
475 :
476 : bool
477 5425 : DisplacedProblem::hasVariable(const std::string & var_name) const
478 : {
479 8780 : for (auto & nl : _displaced_solver_systems)
480 5425 : if (nl->hasVariable(var_name))
481 2070 : return true;
482 3355 : if (_displaced_aux->hasVariable(var_name))
483 3344 : return true;
484 :
485 11 : return false;
486 : }
487 :
488 : const MooseVariableFieldBase &
489 34980 : 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 104940 : return getVariableHelper(tid,
495 : var_name,
496 : expected_var_type,
497 : expected_var_field_type,
498 34980 : _displaced_solver_systems,
499 69960 : *_displaced_aux);
500 : }
501 :
502 : MooseVariable &
503 46468 : DisplacedProblem::getStandardVariable(const THREAD_ID tid, const std::string & var_name)
504 : {
505 92782 : for (auto & nl : _displaced_solver_systems)
506 46468 : if (nl->hasVariable(var_name))
507 154 : return nl->getFieldVariable<Real>(tid, var_name);
508 46314 : if (_displaced_aux->hasVariable(var_name))
509 46314 : 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 2729 : DisplacedProblem::hasScalarVariable(const std::string & var_name) const
552 : {
553 5436 : for (auto & nl : _displaced_solver_systems)
554 2729 : if (nl->hasScalarVariable(var_name))
555 22 : return true;
556 2707 : if (_displaced_aux->hasScalarVariable(var_name))
557 0 : return true;
558 :
559 2707 : return false;
560 : }
561 :
562 : MooseVariableScalar &
563 22 : DisplacedProblem::getScalarVariable(const THREAD_ID tid, const std::string & var_name)
564 : {
565 22 : for (auto & nl : _displaced_solver_systems)
566 22 : if (nl->hasScalarVariable(var_name))
567 22 : 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 13 : DisplacedProblem::getSystem(const std::string & var_name)
576 : {
577 26 : for (const auto sys_num : make_range(_eq.n_systems()))
578 : {
579 26 : auto & sys = _eq.get_system(sys_num);
580 26 : if (sys.has_variable(var_name))
581 13 : return sys;
582 : }
583 :
584 0 : mooseError("Unable to find a system containing the variable " + var_name);
585 : }
586 :
587 : void
588 3813 : DisplacedProblem::addVariable(const std::string & var_type,
589 : const std::string & name,
590 : InputParameters & parameters,
591 : const unsigned int nl_system_number)
592 : {
593 3813 : _displaced_solver_systems[nl_system_number]->addVariable(var_type, name, parameters);
594 3813 : }
595 :
596 : void
597 11073 : DisplacedProblem::addAuxVariable(const std::string & var_type,
598 : const std::string & name,
599 : InputParameters & parameters)
600 : {
601 11073 : _displaced_aux->addVariable(var_type, name, parameters);
602 11073 : }
603 :
604 : unsigned int
605 13006720 : DisplacedProblem::currentNlSysNum() const
606 : {
607 13006720 : return _mproblem.currentNlSysNum();
608 : }
609 :
610 : unsigned int
611 0 : DisplacedProblem::currentLinearSysNum() const
612 : {
613 0 : return _mproblem.currentLinearSysNum();
614 : }
615 :
616 : void
617 8531238 : DisplacedProblem::prepare(const Elem * elem, const THREAD_ID tid)
618 : {
619 17062448 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
620 : {
621 8531238 : _assembly[tid][nl_sys_num]->reinit(elem);
622 8531210 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
623 : // This method is called outside of residual/Jacobian callbacks during initial condition
624 : // evaluation
625 8531210 : if (!_mproblem.hasJacobian() || !_mproblem.constJacobian())
626 8531210 : _assembly[tid][nl_sys_num]->prepareJacobianBlock();
627 8531210 : _assembly[tid][nl_sys_num]->prepareResidual();
628 : }
629 :
630 8531210 : _displaced_aux->prepare(tid);
631 8531210 : }
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 35278 : DisplacedProblem::setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid)
665 : {
666 35278 : SubdomainID did = elem->subdomain_id();
667 70556 : for (auto & assembly : _assembly[tid])
668 35278 : assembly->setCurrentSubdomainID(did);
669 35278 : }
670 :
671 : void
672 143946 : DisplacedProblem::setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid)
673 : {
674 143946 : SubdomainID did = elem->neighbor_ptr(side)->subdomain_id();
675 287892 : for (auto & assembly : _assembly[tid])
676 143946 : assembly->setCurrentNeighborSubdomainID(did);
677 143946 : }
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 57704 : DisplacedProblem::prepareAssembly(const THREAD_ID tid)
691 : {
692 57704 : _assembly[tid][currentNlSysNum()]->prepare();
693 57704 : }
694 :
695 : void
696 57704 : DisplacedProblem::prepareAssemblyNeighbor(const THREAD_ID tid)
697 : {
698 57704 : _assembly[tid][currentNlSysNum()]->prepareNeighbor();
699 57704 : }
700 :
701 : bool
702 5512 : DisplacedProblem::reinitDirac(const Elem * elem, const THREAD_ID tid)
703 : {
704 5512 : std::vector<Point> & points = _dirac_kernel_info.getPoints()[elem].first;
705 :
706 5512 : unsigned int n_points = points.size();
707 :
708 5512 : if (n_points)
709 : {
710 11024 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
711 : {
712 5512 : _assembly[tid][nl_sys_num]->reinitAtPhysical(elem, points);
713 5512 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
714 : }
715 5512 : _displaced_aux->prepare(tid);
716 :
717 5512 : reinitElem(elem, tid);
718 : }
719 :
720 5512 : _assembly[tid][currentNlSysNum()]->prepare();
721 :
722 5512 : return n_points > 0;
723 : }
724 :
725 : void
726 4961974 : DisplacedProblem::reinitElem(const Elem * elem, const THREAD_ID tid)
727 : {
728 9923948 : for (auto & nl : _displaced_solver_systems)
729 4961974 : nl->reinitElem(elem, tid);
730 4961974 : _displaced_aux->reinitElem(elem, tid);
731 4961974 : }
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 94464 : DisplacedProblem::reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid)
754 : {
755 188928 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
756 : {
757 94464 : _assembly[tid][nl_sys_num]->reinit(elem, side);
758 94464 : _displaced_solver_systems[nl_sys_num]->reinitElemFace(elem, side, tid);
759 : }
760 94464 : _displaced_aux->reinitElemFace(elem, side, tid);
761 94464 : }
762 :
763 : void
764 1122418 : DisplacedProblem::reinitNode(const Node * node, const THREAD_ID tid)
765 : {
766 2244836 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
767 : {
768 1122418 : _assembly[tid][nl_sys_num]->reinit(node);
769 1122418 : _displaced_solver_systems[nl_sys_num]->reinitNode(node, tid);
770 : }
771 1122418 : _displaced_aux->reinitNode(node, tid);
772 1122418 : }
773 :
774 : void
775 4893197 : DisplacedProblem::reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid)
776 : {
777 9786394 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
778 : {
779 4893197 : _assembly[tid][nl_sys_num]->reinit(node);
780 4893197 : _displaced_solver_systems[nl_sys_num]->reinitNodeFace(node, bnd_id, tid);
781 : }
782 4893197 : _displaced_aux->reinitNodeFace(node, bnd_id, tid);
783 4893197 : }
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 71940 : DisplacedProblem::reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid)
803 : {
804 71940 : reinitNeighbor(elem, side, tid, nullptr);
805 71940 : }
806 :
807 : void
808 143946 : DisplacedProblem::reinitNeighbor(const Elem * elem,
809 : unsigned int side,
810 : const THREAD_ID tid,
811 : const std::vector<Point> * neighbor_reference_points)
812 : {
813 143946 : setNeighborSubdomainID(elem, side, tid);
814 :
815 143946 : const Elem * neighbor = elem->neighbor_ptr(side);
816 143946 : unsigned int neighbor_side = neighbor->which_neighbor_am_i(elem);
817 :
818 287892 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
819 : {
820 143946 : _assembly[tid][nl_sys_num]->reinitElemAndNeighbor(
821 : elem, side, neighbor, neighbor_side, neighbor_reference_points);
822 143946 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
823 : // Called during stateful material property evaluation outside of solve
824 143946 : _assembly[tid][nl_sys_num]->prepareNeighbor();
825 : }
826 143946 : _displaced_aux->prepareNeighbor(tid);
827 :
828 287892 : for (auto & nl : _displaced_solver_systems)
829 : {
830 143946 : nl->reinitElemFace(elem, side, tid);
831 143946 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
832 : }
833 143946 : _displaced_aux->reinitElemFace(elem, side, tid);
834 143946 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
835 143946 : }
836 :
837 : void
838 57704 : 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 115408 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
847 : {
848 : // Reinit shape functions
849 57704 : _assembly[tid][nl_sys_num]->reinitNeighborAtPhysical(neighbor, neighbor_side, physical_points);
850 :
851 : // Set the neighbor dof indices
852 57704 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
853 : }
854 57704 : _displaced_aux->prepareNeighbor(tid);
855 :
856 57704 : prepareAssemblyNeighbor(tid);
857 :
858 : // Compute values at the points
859 115408 : for (auto & nl : _displaced_solver_systems)
860 57704 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
861 57704 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
862 57704 : }
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 71940 : DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem,
892 : unsigned int side,
893 : const THREAD_ID tid)
894 : {
895 71940 : reinitNeighbor(elem, side, tid);
896 :
897 71940 : const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
898 71940 : if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0)
899 660 : reinitLowerDElem(lower_d_elem, tid);
900 : else
901 : {
902 : // with mesh refinement, lower-dimensional element might be defined on neighbor side
903 71280 : auto & neighbor = _assembly[tid][currentNlSysNum()]->neighbor();
904 71280 : auto & neighbor_side = _assembly[tid][currentNlSysNum()]->neighborSide();
905 71280 : const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side);
906 71280 : if (lower_d_elem_neighbor &&
907 71280 : _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 71940 : }
917 :
918 : void
919 124381 : DisplacedProblem::reinitScalars(const THREAD_ID tid,
920 : bool reinit_for_derivative_reordering /*=false*/)
921 : {
922 248762 : for (auto & nl : _displaced_solver_systems)
923 124381 : nl->reinitScalars(tid, reinit_for_derivative_reordering);
924 124381 : _displaced_aux->reinitScalars(tid, reinit_for_derivative_reordering);
925 124381 : }
926 :
927 : void
928 90 : DisplacedProblem::reinitOffDiagScalars(const THREAD_ID tid)
929 : {
930 90 : _assembly[tid][currentNlSysNum()]->prepareOffDiagScalar();
931 90 : }
932 :
933 : void
934 2539 : DisplacedProblem::getDiracElements(std::set<const Elem *> & elems)
935 : {
936 2539 : elems = _dirac_kernel_info.getElements();
937 2539 : }
938 :
939 : void
940 159166 : DisplacedProblem::clearDiracInfo()
941 : {
942 159166 : _dirac_kernel_info.clearPoints();
943 159166 : }
944 :
945 : void
946 5288 : DisplacedProblem::addResidual(const THREAD_ID tid)
947 : {
948 10576 : _assembly[tid][currentNlSysNum()]->addResidual(Assembly::GlobalDataKey{},
949 5288 : currentResidualVectorTags());
950 5288 : }
951 :
952 : void
953 68460 : DisplacedProblem::addResidualNeighbor(const THREAD_ID tid)
954 : {
955 136920 : _assembly[tid][currentNlSysNum()]->addResidualNeighbor(Assembly::GlobalDataKey{},
956 68460 : currentResidualVectorTags());
957 68460 : }
958 :
959 : void
960 68688 : DisplacedProblem::addResidualLower(const THREAD_ID tid)
961 : {
962 137376 : _assembly[tid][currentNlSysNum()]->addResidualLower(Assembly::GlobalDataKey{},
963 68688 : currentResidualVectorTags());
964 68688 : }
965 :
966 : void
967 48 : DisplacedProblem::addCachedResidualDirectly(NumericVector<Number> & residual, const THREAD_ID tid)
968 : {
969 96 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
970 48 : _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 96 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
977 48 : _displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag()))
978 144 : _assembly[tid][currentNlSysNum()]->addCachedResidualDirectly(
979 : residual,
980 96 : Assembly::GlobalDataKey{},
981 48 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag()));
982 :
983 : // We do this because by adding the cached residual directly, we cannot ensure that all of the
984 : // cached residuals are emptied after only the two add calls above
985 48 : _assembly[tid][currentNlSysNum()]->clearCachedResiduals(Assembly::GlobalDataKey{});
986 48 : }
987 :
988 : void
989 0 : DisplacedProblem::setResidual(NumericVector<Number> & residual, const THREAD_ID tid)
990 : {
991 0 : _assembly[tid][currentNlSysNum()]->setResidual(
992 : residual,
993 0 : Assembly::GlobalDataKey{},
994 0 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->residualVectorTag()));
995 0 : }
996 :
997 : void
998 0 : DisplacedProblem::setResidualNeighbor(NumericVector<Number> & residual, const THREAD_ID tid)
999 : {
1000 0 : _assembly[tid][currentNlSysNum()]->setResidualNeighbor(
1001 : residual,
1002 0 : Assembly::GlobalDataKey{},
1003 0 : getVectorTag(_displaced_solver_systems[currentNlSysNum()]->residualVectorTag()));
1004 0 : }
1005 :
1006 : void
1007 224 : DisplacedProblem::addJacobian(const THREAD_ID tid)
1008 : {
1009 224 : _assembly[tid][currentNlSysNum()]->addJacobian(Assembly::GlobalDataKey{});
1010 224 : }
1011 :
1012 : void
1013 0 : DisplacedProblem::addJacobianNonlocal(const THREAD_ID tid)
1014 : {
1015 0 : _assembly[tid][currentNlSysNum()]->addJacobianNonlocal(Assembly::GlobalDataKey{});
1016 0 : }
1017 :
1018 : void
1019 54 : DisplacedProblem::addJacobianNeighbor(const THREAD_ID tid)
1020 : {
1021 54 : _assembly[tid][currentNlSysNum()]->addJacobianNeighbor(Assembly::GlobalDataKey{});
1022 54 : }
1023 :
1024 : void
1025 3456 : DisplacedProblem::addJacobianNeighborLowerD(const THREAD_ID tid)
1026 : {
1027 3456 : _assembly[tid][currentNlSysNum()]->addJacobianNeighborLowerD(Assembly::GlobalDataKey{});
1028 3456 : }
1029 :
1030 : void
1031 216 : DisplacedProblem::addJacobianLowerD(const THREAD_ID tid)
1032 : {
1033 216 : _assembly[tid][currentNlSysNum()]->addJacobianLowerD(Assembly::GlobalDataKey{});
1034 216 : }
1035 :
1036 : void
1037 0 : DisplacedProblem::cacheJacobianNonlocal(const THREAD_ID tid)
1038 : {
1039 0 : _assembly[tid][currentNlSysNum()]->cacheJacobianNonlocal(Assembly::GlobalDataKey{});
1040 0 : }
1041 :
1042 : void
1043 0 : DisplacedProblem::addJacobianBlockTags(SparseMatrix<Number> & jacobian,
1044 : unsigned int ivar,
1045 : unsigned int jvar,
1046 : const DofMap & dof_map,
1047 : std::vector<dof_id_type> & dof_indices,
1048 : const std::set<TagID> & tags,
1049 : const THREAD_ID tid)
1050 : {
1051 0 : _assembly[tid][currentNlSysNum()]->addJacobianBlockTags(
1052 0 : jacobian, ivar, jvar, dof_map, dof_indices, Assembly::GlobalDataKey{}, tags);
1053 0 : }
1054 :
1055 : void
1056 0 : DisplacedProblem::addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
1057 : unsigned int ivar,
1058 : unsigned int jvar,
1059 : const DofMap & dof_map,
1060 : const std::vector<dof_id_type> & idof_indices,
1061 : const std::vector<dof_id_type> & jdof_indices,
1062 : const std::set<TagID> & tags,
1063 : const THREAD_ID tid)
1064 : {
1065 0 : _assembly[tid][currentNlSysNum()]->addJacobianBlockNonlocalTags(
1066 0 : jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices, Assembly::GlobalDataKey{}, tags);
1067 0 : }
1068 :
1069 : void
1070 0 : DisplacedProblem::addJacobianNeighbor(SparseMatrix<Number> & jacobian,
1071 : unsigned int ivar,
1072 : unsigned int jvar,
1073 : const DofMap & dof_map,
1074 : std::vector<dof_id_type> & dof_indices,
1075 : std::vector<dof_id_type> & neighbor_dof_indices,
1076 : const std::set<TagID> & tags,
1077 : const THREAD_ID tid)
1078 : {
1079 0 : _assembly[tid][currentNlSysNum()]->addJacobianNeighborTags(jacobian,
1080 : ivar,
1081 : jvar,
1082 : dof_map,
1083 : dof_indices,
1084 : neighbor_dof_indices,
1085 0 : Assembly::GlobalDataKey{},
1086 : tags);
1087 0 : }
1088 :
1089 : void
1090 842888 : DisplacedProblem::prepareShapes(unsigned int var, const THREAD_ID tid)
1091 : {
1092 842888 : _assembly[tid][currentNlSysNum()]->copyShapes(var);
1093 842888 : }
1094 :
1095 : void
1096 10056 : DisplacedProblem::prepareFaceShapes(unsigned int var, const THREAD_ID tid)
1097 : {
1098 10056 : _assembly[tid][currentNlSysNum()]->copyFaceShapes(var);
1099 10056 : }
1100 :
1101 : void
1102 3768 : DisplacedProblem::prepareNeighborShapes(unsigned int var, const THREAD_ID tid)
1103 : {
1104 3768 : _assembly[tid][currentNlSysNum()]->copyNeighborShapes(var);
1105 3768 : }
1106 :
1107 : void
1108 4478 : DisplacedProblem::updateGeomSearch(GeometricSearchData::GeometricSearchType type)
1109 : {
1110 4478 : TIME_SECTION("updateGeometricSearch", 3, "Updating Displaced GeometricSearch");
1111 :
1112 4478 : _geometric_search_data.update(type);
1113 4478 : }
1114 :
1115 : void
1116 557 : DisplacedProblem::meshChanged(const bool contract_mesh, const bool clean_refinement_flags)
1117 : {
1118 : // The mesh changed. The displaced equations system object only holds Systems, so calling
1119 : // EquationSystems::reinit only prolongs/restricts the solution vectors, which is something that
1120 : // needs to happen for every step of mesh adaptivity.
1121 557 : _eq.reinit();
1122 557 : if (contract_mesh)
1123 : // Once vectors are restricted, we can delete children of coarsened elements
1124 497 : _mesh.getMesh().contract();
1125 557 : if (clean_refinement_flags)
1126 : {
1127 : // Finally clean refinement flags so that if someone tries to project vectors again without
1128 : // an intervening mesh refinement to clean flags they won't run into trouble
1129 497 : MeshRefinement refinement(_mesh.getMesh());
1130 497 : refinement.clean_refinement_flags();
1131 497 : }
1132 :
1133 : // Since the mesh has changed, we need to make sure that we update any of our
1134 : // MOOSE-system specific data.
1135 1114 : for (auto & nl : _displaced_solver_systems)
1136 557 : nl->reinit();
1137 557 : _displaced_aux->reinit();
1138 :
1139 : // We've performed some mesh adaptivity. We need to
1140 : // clear any quadrature nodes such that when we build the boundary node lists in
1141 : // MooseMesh::meshChanged we don't have any extraneous extra boundary nodes lying around
1142 557 : _mesh.clearQuadratureNodes();
1143 :
1144 557 : _mesh.meshChanged();
1145 :
1146 : // Before performing mesh adaptivity we un-displaced the mesh. We need to re-displace the mesh and
1147 : // then reinitialize GeometricSearchData such that we have all the correct geometric information
1148 : // for the changed mesh
1149 557 : updateMesh(/*mesh_changing=*/true);
1150 557 : }
1151 :
1152 : void
1153 893059 : DisplacedProblem::addGhostedElem(dof_id_type elem_id)
1154 : {
1155 893059 : _mproblem.addGhostedElem(elem_id);
1156 893059 : }
1157 :
1158 : void
1159 28798 : DisplacedProblem::addGhostedBoundary(BoundaryID boundary_id)
1160 : {
1161 28798 : _mproblem.addGhostedBoundary(boundary_id);
1162 28798 : }
1163 :
1164 : void
1165 0 : DisplacedProblem::ghostGhostedBoundaries()
1166 : {
1167 0 : _mproblem.ghostGhostedBoundaries();
1168 0 : }
1169 :
1170 : MooseMesh &
1171 454525 : DisplacedProblem::refMesh()
1172 : {
1173 454525 : return _ref_mesh;
1174 : }
1175 :
1176 : bool
1177 0 : DisplacedProblem::solverSystemConverged(const unsigned int sys_num)
1178 : {
1179 0 : return _mproblem.converged(sys_num);
1180 : }
1181 :
1182 : bool
1183 0 : DisplacedProblem::computingPreSMOResidual(const unsigned int nl_sys_num) const
1184 : {
1185 0 : return _mproblem.computingPreSMOResidual(nl_sys_num);
1186 : }
1187 :
1188 : void
1189 0 : DisplacedProblem::onTimestepBegin()
1190 : {
1191 0 : }
1192 :
1193 : void
1194 0 : DisplacedProblem::onTimestepEnd()
1195 : {
1196 0 : }
1197 :
1198 : void
1199 633 : DisplacedProblem::undisplaceMesh()
1200 : {
1201 : // If undisplaceMesh() is called during initial adaptivity, it is
1202 : // not valid to call _mesh.getActiveSemiLocalNodeRange() since it is
1203 : // not set up yet. So we are creating the Range by hand.
1204 : //
1205 : // We must undisplace *all* our nodes to the _ref_mesh
1206 : // configuration, not just the local ones, since the partitioners
1207 : // require this. We are using the GRAIN_SIZE=1 from MooseMesh.C,
1208 : // not sure how this value was decided upon.
1209 : //
1210 : // (DRG: The grainsize parameter is ultimately passed to TBB to help
1211 : // it choose how to split up the range. A grainsize of 1 says "split
1212 : // it as much as you want". Years ago I experimentally found that it
1213 : // didn't matter much and that using 1 was fine.)
1214 : //
1215 : // Note: we don't have to invalidate/update as much stuff as
1216 : // DisplacedProblem::updateMesh() does, since this will be handled
1217 : // by a later call to updateMesh().
1218 1266 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
1219 1266 : _mesh.getMesh().nodes_end(),
1220 633 : /*grainsize=*/1);
1221 :
1222 633 : ResetDisplacedMeshThread rdmt(_mproblem, *this);
1223 :
1224 : // Undisplace the mesh using threads.
1225 633 : Threads::parallel_reduce(node_range, rdmt);
1226 633 : }
1227 :
1228 : LineSearch *
1229 0 : DisplacedProblem::getLineSearch()
1230 : {
1231 0 : return _mproblem.getLineSearch();
1232 : }
1233 :
1234 : const CouplingMatrix *
1235 0 : DisplacedProblem::couplingMatrix(const unsigned int nl_sys_num) const
1236 : {
1237 0 : return _mproblem.couplingMatrix(nl_sys_num);
1238 : }
1239 :
1240 : bool
1241 566617 : DisplacedProblem::computingScalingJacobian() const
1242 : {
1243 566617 : return _mproblem.computingScalingJacobian();
1244 : }
1245 :
1246 : bool
1247 0 : DisplacedProblem::computingScalingResidual() const
1248 : {
1249 0 : return _mproblem.computingScalingResidual();
1250 : }
1251 :
1252 : void
1253 2221 : DisplacedProblem::initialSetup()
1254 : {
1255 2221 : SubProblem::initialSetup();
1256 :
1257 4442 : for (auto & nl : _displaced_solver_systems)
1258 2221 : nl->initialSetup();
1259 2221 : _displaced_aux->initialSetup();
1260 2221 : }
1261 :
1262 : void
1263 33478 : DisplacedProblem::timestepSetup()
1264 : {
1265 33478 : SubProblem::timestepSetup();
1266 :
1267 66956 : for (auto & nl : _displaced_solver_systems)
1268 33478 : nl->timestepSetup();
1269 33478 : _displaced_aux->timestepSetup();
1270 33478 : }
1271 :
1272 : void
1273 158722 : DisplacedProblem::customSetup(const ExecFlagType & exec_type)
1274 : {
1275 158722 : SubProblem::customSetup(exec_type);
1276 :
1277 317444 : for (auto & nl : _displaced_solver_systems)
1278 158722 : nl->customSetup(exec_type);
1279 158722 : _displaced_aux->customSetup(exec_type);
1280 158722 : }
1281 :
1282 : void
1283 137122 : DisplacedProblem::residualSetup()
1284 : {
1285 137122 : SubProblem::residualSetup();
1286 :
1287 274244 : for (auto & nl : _displaced_solver_systems)
1288 137122 : nl->residualSetup();
1289 137122 : _displaced_aux->residualSetup();
1290 137122 : }
1291 :
1292 : void
1293 23319 : DisplacedProblem::jacobianSetup()
1294 : {
1295 23319 : SubProblem::jacobianSetup();
1296 :
1297 46638 : for (auto & nl : _displaced_solver_systems)
1298 23319 : nl->jacobianSetup();
1299 23319 : _displaced_aux->jacobianSetup();
1300 23319 : }
1301 :
1302 : void
1303 300 : DisplacedProblem::haveADObjects(const bool have_ad_objects)
1304 : {
1305 300 : _have_ad_objects = have_ad_objects;
1306 300 : _mproblem.SubProblem::haveADObjects(have_ad_objects);
1307 300 : }
1308 :
1309 : std::pair<bool, unsigned int>
1310 34980 : DisplacedProblem::determineSolverSystem(const std::string & var_name,
1311 : const bool error_if_not_found) const
1312 : {
1313 34980 : return _mproblem.determineSolverSystem(var_name, error_if_not_found);
1314 : }
1315 :
1316 : Assembly &
1317 59534016 : DisplacedProblem::assembly(const THREAD_ID tid, const unsigned int sys_num)
1318 : {
1319 : mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
1320 : mooseAssert(sys_num < _assembly[tid].size(),
1321 : "System number larger than the assembly container size");
1322 59534016 : return *_assembly[tid][sys_num];
1323 : }
1324 :
1325 : const Assembly &
1326 48604 : DisplacedProblem::assembly(const THREAD_ID tid, const unsigned int sys_num) const
1327 : {
1328 : mooseAssert(tid < _assembly.size(), "Assembly objects not initialized");
1329 : mooseAssert(sys_num < _assembly[tid].size(),
1330 : "System number larger than the assembly container size");
1331 48604 : return *_assembly[tid][sys_num];
1332 : }
1333 :
1334 : std::size_t
1335 10354426 : DisplacedProblem::numNonlinearSystems() const
1336 : {
1337 10354426 : return _mproblem.numNonlinearSystems();
1338 : }
1339 :
1340 : std::size_t
1341 0 : DisplacedProblem::numLinearSystems() const
1342 : {
1343 0 : return _mproblem.numLinearSystems();
1344 : }
1345 :
1346 : std::size_t
1347 398080 : DisplacedProblem::numSolverSystems() const
1348 : {
1349 398080 : return _mproblem.numSolverSystems();
1350 : }
1351 :
1352 : const std::vector<VectorTag> &
1353 9342862 : DisplacedProblem::currentResidualVectorTags() const
1354 : {
1355 9342862 : return _mproblem.currentResidualVectorTags();
1356 : }
1357 :
1358 : bool
1359 61840539 : DisplacedProblem::safeAccessTaggedMatrices() const
1360 : {
1361 61840539 : return _mproblem.safeAccessTaggedMatrices();
1362 : }
1363 :
1364 : bool
1365 490 : DisplacedProblem::safeAccessTaggedVectors() const
1366 : {
1367 490 : return _mproblem.safeAccessTaggedVectors();
1368 : }
1369 :
1370 : void
1371 0 : DisplacedProblem::needFV()
1372 : {
1373 0 : _mproblem.needFV();
1374 0 : }
1375 :
1376 : bool
1377 326239 : DisplacedProblem::haveFV() const
1378 : {
1379 326239 : return _mproblem.haveFV();
1380 : }
1381 :
1382 : bool
1383 1766776 : DisplacedProblem::hasNonlocalCoupling() const
1384 : {
1385 1766776 : return _mproblem.hasNonlocalCoupling();
1386 : }
1387 :
1388 : unsigned int
1389 0 : DisplacedProblem::nlSysNum(const NonlinearSystemName & nl_sys_name) const
1390 : {
1391 0 : return _mproblem.nlSysNum(nl_sys_name);
1392 : }
1393 :
1394 : unsigned int
1395 0 : DisplacedProblem::linearSysNum(const LinearSystemName & sys_name) const
1396 : {
1397 0 : return _mproblem.linearSysNum(sys_name);
1398 : }
1399 :
1400 : unsigned int
1401 0 : DisplacedProblem::solverSysNum(const SolverSystemName & sys_name) const
1402 : {
1403 0 : return _mproblem.solverSysNum(sys_name);
1404 : }
1405 :
1406 : const libMesh::CouplingMatrix &
1407 2432 : DisplacedProblem::nonlocalCouplingMatrix(const unsigned i) const
1408 : {
1409 2432 : return _mproblem.nonlocalCouplingMatrix(i);
1410 : }
1411 :
1412 : bool
1413 0 : DisplacedProblem::checkNonlocalCouplingRequirement() const
1414 : {
1415 0 : return _mproblem.checkNonlocalCouplingRequirement();
1416 : }
1417 :
1418 163565 : DisplacedProblem::UpdateDisplacedMeshThread::UpdateDisplacedMeshThread(
1419 163565 : FEProblemBase & fe_problem, DisplacedProblem & displaced_problem)
1420 : : ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>(fe_problem),
1421 163565 : _displaced_problem(displaced_problem),
1422 327130 : _ref_mesh(_displaced_problem.refMesh()),
1423 163565 : _nl_soln(_displaced_problem._nl_solution),
1424 163565 : _aux_soln(*_displaced_problem._aux_solution),
1425 163565 : _has_displacement(false)
1426 : {
1427 163565 : this->init();
1428 163565 : }
1429 :
1430 14492 : DisplacedProblem::UpdateDisplacedMeshThread::UpdateDisplacedMeshThread(
1431 14492 : UpdateDisplacedMeshThread & x, Threads::split split)
1432 : : ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>(x, split),
1433 14492 : _displaced_problem(x._displaced_problem),
1434 14492 : _ref_mesh(x._ref_mesh),
1435 14492 : _nl_soln(x._nl_soln),
1436 14492 : _aux_soln(x._aux_soln),
1437 14492 : _sys_to_nonghost_and_ghost_soln(x._sys_to_nonghost_and_ghost_soln),
1438 14492 : _sys_to_var_num_and_direction(x._sys_to_var_num_and_direction),
1439 14492 : _has_displacement(x._has_displacement)
1440 : {
1441 14492 : }
1442 :
1443 : void
1444 163565 : DisplacedProblem::UpdateDisplacedMeshThread::init()
1445 : {
1446 163565 : std::vector<std::string> & displacement_variables = _displaced_problem._displacements;
1447 163565 : unsigned int num_displacements = displacement_variables.size();
1448 163565 : auto & es = _displaced_problem.es();
1449 :
1450 163565 : _sys_to_var_num_and_direction.clear();
1451 163565 : _sys_to_nonghost_and_ghost_soln.clear();
1452 :
1453 546260 : for (unsigned int i = 0; i < num_displacements; i++)
1454 : {
1455 382695 : std::string displacement_name = displacement_variables[i];
1456 :
1457 501881 : for (const auto sys_num : make_range(es.n_systems()))
1458 : {
1459 501881 : auto & sys = es.get_system(sys_num);
1460 501881 : if (sys.has_variable(displacement_name))
1461 : {
1462 382695 : auto & val = _sys_to_var_num_and_direction[sys.number()];
1463 382695 : val.first.push_back(sys.variable_number(displacement_name));
1464 382695 : val.second.push_back(i);
1465 382695 : break;
1466 : }
1467 : }
1468 382695 : }
1469 :
1470 329750 : for (const auto & pr : _sys_to_var_num_and_direction)
1471 : {
1472 166185 : auto & sys = es.get_system(pr.first);
1473 : mooseAssert(sys.number() <= _nl_soln.size(),
1474 : "The system number should always be less than or equal to the number of nonlinear "
1475 : "systems. If it is equal, then this system is the auxiliary system");
1476 : const NumericVector<Number> * const nonghost_soln =
1477 166185 : sys.number() < _nl_soln.size() ? _nl_soln[sys.number()] : &_aux_soln;
1478 332370 : _sys_to_nonghost_and_ghost_soln.emplace(
1479 166185 : sys.number(),
1480 0 : std::make_pair(nonghost_soln,
1481 332370 : NumericVector<Number>::build(nonghost_soln->comm()).release()));
1482 : }
1483 :
1484 163565 : ConstNodeRange node_range(_ref_mesh.getMesh().nodes_begin(), _ref_mesh.getMesh().nodes_end());
1485 :
1486 329750 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1487 : {
1488 166185 : auto & sys = es.get_system(sys_num);
1489 : AllNodesSendListThread send_list(
1490 166185 : this->_fe_problem, _ref_mesh, var_num_and_direction.first, sys);
1491 166185 : Threads::parallel_reduce(node_range, send_list);
1492 166185 : send_list.unique();
1493 166185 : auto & [soln, ghost_soln] = libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num);
1494 332370 : ghost_soln->init(
1495 166185 : soln->size(), soln->local_size(), send_list.send_list(), true, libMesh::GHOSTED);
1496 166185 : soln->localize(*ghost_soln, send_list.send_list());
1497 166185 : }
1498 :
1499 163565 : _has_displacement = false;
1500 163565 : }
1501 :
1502 : void
1503 19967808 : DisplacedProblem::UpdateDisplacedMeshThread::onNode(NodeRange::const_iterator & nd)
1504 : {
1505 19967808 : Node & displaced_node = *(*nd);
1506 :
1507 19967808 : Node & reference_node = _ref_mesh.nodeRef(displaced_node.id());
1508 :
1509 40474144 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1510 : {
1511 20506336 : auto & var_numbers = var_num_and_direction.first;
1512 20506336 : auto & directions = var_num_and_direction.second;
1513 62972930 : for (const auto i : index_range(var_numbers))
1514 : {
1515 42466594 : const auto direction = directions[i];
1516 42466594 : if (reference_node.n_dofs(sys_num, var_numbers[i]) > 0)
1517 : {
1518 42466454 : Real coord = reference_node(direction) +
1519 84932908 : (*libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num).second)(
1520 42466454 : reference_node.dof_number(sys_num, var_numbers[i], 0));
1521 42466454 : if (displaced_node(direction) != coord)
1522 : {
1523 3617366 : displaced_node(direction) = coord;
1524 3617366 : _has_displacement = true;
1525 : }
1526 : }
1527 : }
1528 : }
1529 19967808 : }
|
