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 18391 : DisplacedProblem::validParams()
35 : {
36 18391 : InputParameters params = SubProblem::validParams();
37 18391 : params.addClassDescription(
38 : "A Problem object for providing access to the displaced finite element "
39 : "mesh and associated variables.");
40 18391 : params.addPrivateParam<MooseMesh *>("mesh");
41 18391 : params.addPrivateParam<std::vector<std::string>>("displacements", {});
42 18391 : return params;
43 0 : }
44 :
45 2063 : DisplacedProblem::DisplacedProblem(const InputParameters & parameters)
46 : : SubProblem(parameters),
47 2063 : _mproblem(parameters.have_parameter<FEProblemBase *>("_fe_problem_base")
48 4126 : ? *getParam<FEProblemBase *>("_fe_problem_base")
49 2063 : : *getParam<FEProblem *>("_fe_problem")),
50 2063 : _mesh(*getParam<MooseMesh *>("mesh")),
51 2063 : _eq(_mesh),
52 2063 : _ref_mesh(_mproblem.mesh()),
53 2063 : _displacements(getParam<std::vector<std::string>>("displacements")),
54 4126 : _geometric_search_data(*this, _mesh)
55 :
56 : {
57 : // Disable refinement/coarsening in EquationSystems::reinit because we already do this ourselves
58 2063 : _eq.disable_refine_in_reinit();
59 :
60 : // TODO: Move newAssemblyArray further up to SubProblem so that we can use it here
61 2063 : unsigned int n_threads = libMesh::n_threads();
62 :
63 2063 : _assembly.resize(n_threads);
64 4126 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
65 : {
66 2063 : _displaced_solver_systems.emplace_back(std::make_unique<DisplacedSystem>(
67 : *this,
68 : _mproblem,
69 2063 : _mproblem.getNonlinearSystemBase(nl_sys_num),
70 4126 : "displaced_" + _mproblem.getNonlinearSystemBase(nl_sys_num).name() + "_" +
71 2063 : std::to_string(nl_sys_num),
72 2063 : Moose::VAR_SOLVER));
73 2063 : auto & displaced_nl = _displaced_solver_systems.back();
74 :
75 4325 : for (unsigned int i = 0; i < n_threads; ++i)
76 2262 : _assembly[i].emplace_back(std::make_unique<Assembly>(*displaced_nl, i));
77 : }
78 :
79 2063 : _nl_solution.resize(_displaced_solver_systems.size(), nullptr);
80 :
81 : _displaced_aux =
82 4126 : std::make_unique<DisplacedSystem>(*this,
83 : _mproblem,
84 2063 : _mproblem.getAuxiliarySystem(),
85 4126 : "displaced_" + _mproblem.getAuxiliarySystem().name(),
86 4126 : 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 2063 : automaticScaling(_mproblem.automaticScaling());
100 :
101 2063 : _mesh.setCoordData(_ref_mesh);
102 2063 : }
103 :
104 4046 : DisplacedProblem::~DisplacedProblem() = default;
105 :
106 : bool
107 110782424 : DisplacedProblem::isTransient() const
108 : {
109 110782424 : return _mproblem.isTransient();
110 : }
111 :
112 : std::set<dof_id_type> &
113 8630 : DisplacedProblem::ghostedElems()
114 : {
115 8630 : return _mproblem.ghostedElems();
116 : }
117 :
118 : void
119 2063 : 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 4325 : for (unsigned int tid = 0; tid < libMesh::n_threads(); ++tid)
127 4524 : for (const auto sys_num : index_range(_assembly[tid]))
128 2262 : _assembly[tid][sys_num]->createQRules(
129 : type, order, volume_order, face_order, block, allow_negative_qweights);
130 2063 : }
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 2063 : DisplacedProblem::init()
150 : {
151 4325 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
152 : {
153 4524 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
154 2262 : _assembly[tid][nl_sys_num]->init(_mproblem.couplingMatrix(nl_sys_num));
155 :
156 4524 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
157 : {
158 2262 : std::vector<std::pair<unsigned int, unsigned short>> disp_numbers_and_directions;
159 7143 : for (const auto direction : index_range(_displacements))
160 : {
161 4881 : const auto & disp_string = _displacements[direction];
162 4881 : const auto & disp_variable = getVariable(tid, disp_string);
163 4881 : if (disp_variable.sys().number() == nl_sys_num)
164 2189 : disp_numbers_and_directions.push_back(
165 4378 : std::make_pair(disp_variable.number(), cast_int<unsigned short>(direction)));
166 : }
167 2262 : _assembly[tid][nl_sys_num]->assignDisplacements(std::move(disp_numbers_and_directions));
168 2262 : }
169 : }
170 :
171 4126 : for (auto & nl : _displaced_solver_systems)
172 : {
173 2063 : nl->dofMap().attach_extra_send_list_function(&extraSendList, nl.get());
174 2063 : nl->preInit();
175 : }
176 :
177 2063 : _displaced_aux->dofMap().attach_extra_send_list_function(&extraSendList, _displaced_aux.get());
178 2063 : _displaced_aux->preInit();
179 :
180 : {
181 2063 : TIME_SECTION("eq::init", 2, "Initializing Displaced Equation System");
182 2063 : _eq.init();
183 2063 : }
184 :
185 4126 : for (auto & nl : _displaced_solver_systems)
186 2063 : nl->postInit();
187 2063 : _displaced_aux->postInit();
188 :
189 2063 : _mesh.meshChanged();
190 :
191 2063 : if (haveFV())
192 12 : _mesh.setupFiniteVolumeMeshData();
193 2063 : }
194 :
195 : void
196 130 : DisplacedProblem::initAdaptivity()
197 : {
198 130 : }
199 :
200 : void
201 1746 : DisplacedProblem::addTimeIntegrator()
202 : {
203 3492 : for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems()))
204 1746 : _displaced_solver_systems[nl_sys_num]->copyTimeIntegrators(
205 1746 : _mproblem.getNonlinearSystemBase(nl_sys_num));
206 1746 : _displaced_aux->copyTimeIntegrators(_mproblem.getAuxiliarySystem());
207 1746 : }
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 489999 : DisplacedProblem::syncAuxSolution(const NumericVector<Number> & aux_soln)
227 : {
228 489999 : (*_displaced_aux->sys().solution) = aux_soln;
229 489999 : _displaced_aux->update();
230 489999 : }
231 :
232 : void
233 310826 : DisplacedProblem::syncSolutions()
234 : {
235 310826 : TIME_SECTION("syncSolutions", 5, "Syncing Displaced Solutions");
236 :
237 621652 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
238 : {
239 310826 : 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 310826 : (*displaced_nl->sys().solution) =
244 621652 : *_mproblem.getNonlinearSystemBase(displaced_nl->number()).currentSolution();
245 310826 : displaced_nl->update();
246 : }
247 310826 : syncAuxSolution(*_mproblem.getAuxiliarySystem().currentSolution());
248 310826 : }
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 150580 : DisplacedProblem::updateMesh(bool mesh_changing)
267 : {
268 150580 : 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 150580 : if (!mesh_changing)
275 150071 : syncSolutions();
276 :
277 301160 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
278 150580 : _nl_solution[nl_sys_num] = _displaced_solver_systems[nl_sys_num]->sys().solution.get();
279 150580 : _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 150580 : if (_mesh.getMesh().is_serial() && !this->refMesh().getMesh().is_serial())
287 0 : this->refMesh().getMesh().allgather();
288 :
289 150580 : if (_mesh.getMesh().is_serial_on_zero() && !this->refMesh().getMesh().is_serial_on_zero())
290 0 : this->refMesh().getMesh().gather_to_zero();
291 :
292 150580 : UpdateDisplacedMeshThread udmt(_mproblem, *this);
293 :
294 : // We displace all nodes, not just semilocal nodes, because
295 : // parallel-inconsistent mesh geometry makes libMesh cry.
296 301160 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
297 301160 : _mesh.getMesh().nodes_end(),
298 150580 : /*grainsize=*/1);
299 :
300 150580 : Threads::parallel_reduce(node_range, udmt);
301 : // Displacement of the mesh has invalidated the point locator data (e.g. bounding boxes)
302 150580 : _mesh.getMesh().clear_point_locator();
303 :
304 : // The mesh has changed. Face information normals, areas, etc. must be re-calculated
305 150580 : 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 150580 : if (mesh_changing)
317 509 : _geometric_search_data.reinit();
318 : else
319 150071 : _geometric_search_data.update();
320 : }
321 0 : catch (MooseException & e)
322 : {
323 0 : _mproblem.setException(e.what());
324 0 : }
325 :
326 150576 : if (udmt.hasDisplacement())
327 64965 : _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 150576 : _mproblem.checkExceptionAndStopSolve(/*print_message=*/false);
333 :
334 : // Since the Mesh changed, update the PointLocator object used by DiracKernels.
335 150576 : _dirac_kernel_info.updatePointLocator(_mesh);
336 150576 : }
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 5749509 : DisplacedProblem::getVectorTag(const TagID tag_id) const
394 : {
395 5749509 : return _mproblem.getVectorTag(tag_id);
396 : }
397 :
398 : TagID
399 107139 : DisplacedProblem::getVectorTagID(const TagName & tag_name) const
400 : {
401 107139 : 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 223 : DisplacedProblem::vectorTagExists(const TagID tag_id) const
412 : {
413 223 : 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 44951 : DisplacedProblem::numVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
424 : {
425 44951 : return _mproblem.numVectorTags(type);
426 : }
427 :
428 : const std::vector<VectorTag> &
429 2262 : DisplacedProblem::getVectorTags(const Moose::VectorTagType type /* = Moose::VECTOR_TAG_ANY */) const
430 : {
431 2262 : return _mproblem.getVectorTags(type);
432 : }
433 :
434 : Moose::VectorTagType
435 113853375 : DisplacedProblem::vectorTagType(const TagID tag_id) const
436 : {
437 113853375 : 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 1509 : DisplacedProblem::getMatrixTagID(const TagName & tag_name) const
448 : {
449 1509 : 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 365 : DisplacedProblem::matrixTagExists(TagID tag_id) const
466 : {
467 365 : return _mproblem.matrixTagExists(tag_id);
468 : }
469 :
470 : unsigned int
471 47213 : DisplacedProblem::numMatrixTags() const
472 : {
473 47213 : return _mproblem.numMatrixTags();
474 : }
475 :
476 : bool
477 5060 : DisplacedProblem::hasVariable(const std::string & var_name) const
478 : {
479 8172 : for (auto & nl : _displaced_solver_systems)
480 5060 : if (nl->hasVariable(var_name))
481 1948 : return true;
482 3112 : if (_displaced_aux->hasVariable(var_name))
483 3102 : return true;
484 :
485 10 : return false;
486 : }
487 :
488 : const MooseVariableFieldBase &
489 32547 : 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 97641 : return getVariableHelper(tid,
495 : var_name,
496 : expected_var_type,
497 : expected_var_field_type,
498 32547 : _displaced_solver_systems,
499 65094 : *_displaced_aux);
500 : }
501 :
502 : MooseVariable &
503 42812 : DisplacedProblem::getStandardVariable(const THREAD_ID tid, const std::string & var_name)
504 : {
505 85480 : for (auto & nl : _displaced_solver_systems)
506 42812 : if (nl->hasVariable(var_name))
507 144 : return nl->getFieldVariable<Real>(tid, var_name);
508 42668 : if (_displaced_aux->hasVariable(var_name))
509 42668 : 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 2545 : DisplacedProblem::hasScalarVariable(const std::string & var_name) const
552 : {
553 5070 : for (auto & nl : _displaced_solver_systems)
554 2545 : if (nl->hasScalarVariable(var_name))
555 20 : return true;
556 2525 : if (_displaced_aux->hasScalarVariable(var_name))
557 0 : return true;
558 :
559 2525 : return false;
560 : }
561 :
562 : MooseVariableScalar &
563 20 : DisplacedProblem::getScalarVariable(const THREAD_ID tid, const std::string & var_name)
564 : {
565 20 : for (auto & nl : _displaced_solver_systems)
566 20 : if (nl->hasScalarVariable(var_name))
567 20 : 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 3532 : DisplacedProblem::addVariable(const std::string & var_type,
589 : const std::string & name,
590 : InputParameters & parameters,
591 : const unsigned int nl_system_number)
592 : {
593 3532 : _displaced_solver_systems[nl_system_number]->addVariable(var_type, name, parameters);
594 3532 : }
595 :
596 : void
597 10206 : DisplacedProblem::addAuxVariable(const std::string & var_type,
598 : const std::string & name,
599 : InputParameters & parameters)
600 : {
601 10206 : _displaced_aux->addVariable(var_type, name, parameters);
602 10206 : }
603 :
604 : unsigned int
605 11587801 : DisplacedProblem::currentNlSysNum() const
606 : {
607 11587801 : return _mproblem.currentNlSysNum();
608 : }
609 :
610 : unsigned int
611 0 : DisplacedProblem::currentLinearSysNum() const
612 : {
613 0 : return _mproblem.currentLinearSysNum();
614 : }
615 :
616 : void
617 7617850 : DisplacedProblem::prepare(const Elem * elem, const THREAD_ID tid)
618 : {
619 15235672 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
620 : {
621 7617850 : _assembly[tid][nl_sys_num]->reinit(elem);
622 7617822 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
623 : // This method is called outside of residual/Jacobian callbacks during initial condition
624 : // evaluation
625 7617822 : if (!_mproblem.hasJacobian() || !_mproblem.constJacobian())
626 7617822 : _assembly[tid][nl_sys_num]->prepareJacobianBlock();
627 7617822 : _assembly[tid][nl_sys_num]->prepareResidual();
628 : }
629 :
630 7617822 : _displaced_aux->prepare(tid);
631 7617822 : }
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 31596 : DisplacedProblem::setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid)
665 : {
666 31596 : SubdomainID did = elem->subdomain_id();
667 63192 : for (auto & assembly : _assembly[tid])
668 31596 : assembly->setCurrentSubdomainID(did);
669 31596 : }
670 :
671 : void
672 129756 : DisplacedProblem::setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid)
673 : {
674 129756 : SubdomainID did = elem->neighbor_ptr(side)->subdomain_id();
675 259512 : for (auto & assembly : _assembly[tid])
676 129756 : assembly->setCurrentNeighborSubdomainID(did);
677 129756 : }
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 51264 : DisplacedProblem::prepareAssembly(const THREAD_ID tid)
691 : {
692 51264 : _assembly[tid][currentNlSysNum()]->prepare();
693 51264 : }
694 :
695 : void
696 51264 : DisplacedProblem::prepareAssemblyNeighbor(const THREAD_ID tid)
697 : {
698 51264 : _assembly[tid][currentNlSysNum()]->prepareNeighbor();
699 51264 : }
700 :
701 : bool
702 4804 : DisplacedProblem::reinitDirac(const Elem * elem, const THREAD_ID tid)
703 : {
704 4804 : std::vector<Point> & points = _dirac_kernel_info.getPoints()[elem].first;
705 :
706 4804 : unsigned int n_points = points.size();
707 :
708 4804 : if (n_points)
709 : {
710 9608 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
711 : {
712 4804 : _assembly[tid][nl_sys_num]->reinitAtPhysical(elem, points);
713 4804 : _displaced_solver_systems[nl_sys_num]->prepare(tid);
714 : }
715 4804 : _displaced_aux->prepare(tid);
716 :
717 4804 : reinitElem(elem, tid);
718 : }
719 :
720 4804 : _assembly[tid][currentNlSysNum()]->prepare();
721 :
722 4804 : return n_points > 0;
723 : }
724 :
725 : void
726 4422879 : DisplacedProblem::reinitElem(const Elem * elem, const THREAD_ID tid)
727 : {
728 8845758 : for (auto & nl : _displaced_solver_systems)
729 4422879 : nl->reinitElem(elem, tid);
730 4422879 : _displaced_aux->reinitElem(elem, tid);
731 4422879 : }
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 83886 : DisplacedProblem::reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid)
754 : {
755 167772 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
756 : {
757 83886 : _assembly[tid][nl_sys_num]->reinit(elem, side);
758 83886 : _displaced_solver_systems[nl_sys_num]->reinitElemFace(elem, side, tid);
759 : }
760 83886 : _displaced_aux->reinitElemFace(elem, side, tid);
761 83886 : }
762 :
763 : void
764 1002259 : DisplacedProblem::reinitNode(const Node * node, const THREAD_ID tid)
765 : {
766 2004518 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
767 : {
768 1002259 : _assembly[tid][nl_sys_num]->reinit(node);
769 1002259 : _displaced_solver_systems[nl_sys_num]->reinitNode(node, tid);
770 : }
771 1002259 : _displaced_aux->reinitNode(node, tid);
772 1002259 : }
773 :
774 : void
775 4339042 : DisplacedProblem::reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid)
776 : {
777 8678084 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
778 : {
779 4339042 : _assembly[tid][nl_sys_num]->reinit(node);
780 4339042 : _displaced_solver_systems[nl_sys_num]->reinitNodeFace(node, bnd_id, tid);
781 : }
782 4339042 : _displaced_aux->reinitNodeFace(node, bnd_id, tid);
783 4339042 : }
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 64848 : DisplacedProblem::reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid)
803 : {
804 64848 : reinitNeighbor(elem, side, tid, nullptr);
805 64848 : }
806 :
807 : void
808 129756 : DisplacedProblem::reinitNeighbor(const Elem * elem,
809 : unsigned int side,
810 : const THREAD_ID tid,
811 : const std::vector<Point> * neighbor_reference_points)
812 : {
813 129756 : setNeighborSubdomainID(elem, side, tid);
814 :
815 129756 : const Elem * neighbor = elem->neighbor_ptr(side);
816 129756 : unsigned int neighbor_side = neighbor->which_neighbor_am_i(elem);
817 :
818 259512 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
819 : {
820 129756 : _assembly[tid][nl_sys_num]->reinitElemAndNeighbor(
821 : elem, side, neighbor, neighbor_side, neighbor_reference_points);
822 129756 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
823 : // Called during stateful material property evaluation outside of solve
824 129756 : _assembly[tid][nl_sys_num]->prepareNeighbor();
825 : }
826 129756 : _displaced_aux->prepareNeighbor(tid);
827 :
828 259512 : for (auto & nl : _displaced_solver_systems)
829 : {
830 129756 : nl->reinitElemFace(elem, side, tid);
831 129756 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
832 : }
833 129756 : _displaced_aux->reinitElemFace(elem, side, tid);
834 129756 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
835 129756 : }
836 :
837 : void
838 51264 : 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 102528 : for (const auto nl_sys_num : index_range(_displaced_solver_systems))
847 : {
848 : // Reinit shape functions
849 51264 : _assembly[tid][nl_sys_num]->reinitNeighborAtPhysical(neighbor, neighbor_side, physical_points);
850 :
851 : // Set the neighbor dof indices
852 51264 : _displaced_solver_systems[nl_sys_num]->prepareNeighbor(tid);
853 : }
854 51264 : _displaced_aux->prepareNeighbor(tid);
855 :
856 51264 : prepareAssemblyNeighbor(tid);
857 :
858 : // Compute values at the points
859 102528 : for (auto & nl : _displaced_solver_systems)
860 51264 : nl->reinitNeighborFace(neighbor, neighbor_side, tid);
861 51264 : _displaced_aux->reinitNeighborFace(neighbor, neighbor_side, tid);
862 51264 : }
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 64848 : DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem,
892 : unsigned int side,
893 : const THREAD_ID tid)
894 : {
895 64848 : reinitNeighbor(elem, side, tid);
896 :
897 64848 : const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
898 64848 : if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0)
899 588 : 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 64848 : }
917 :
918 : void
919 118893 : DisplacedProblem::reinitScalars(const THREAD_ID tid,
920 : bool reinit_for_derivative_reordering /*=false*/)
921 : {
922 237786 : for (auto & nl : _displaced_solver_systems)
923 118893 : nl->reinitScalars(tid, reinit_for_derivative_reordering);
924 118893 : _displaced_aux->reinitScalars(tid, reinit_for_derivative_reordering);
925 118893 : }
926 :
927 : void
928 80 : DisplacedProblem::reinitOffDiagScalars(const THREAD_ID tid)
929 : {
930 80 : _assembly[tid][currentNlSysNum()]->prepareOffDiagScalar();
931 80 : }
932 :
933 : void
934 2316 : DisplacedProblem::getDiracElements(std::set<const Elem *> & elems)
935 : {
936 2316 : elems = _dirac_kernel_info.getElements();
937 2316 : }
938 :
939 : void
940 146566 : DisplacedProblem::clearDiracInfo()
941 : {
942 146566 : _dirac_kernel_info.clearPoints();
943 146566 : }
944 :
945 : void
946 4612 : DisplacedProblem::addResidual(const THREAD_ID tid)
947 : {
948 9224 : _assembly[tid][currentNlSysNum()]->addResidual(Assembly::GlobalDataKey{},
949 4612 : currentResidualVectorTags());
950 4612 : }
951 :
952 : void
953 61756 : DisplacedProblem::addResidualNeighbor(const THREAD_ID tid)
954 : {
955 123512 : _assembly[tid][currentNlSysNum()]->addResidualNeighbor(Assembly::GlobalDataKey{},
956 61756 : currentResidualVectorTags());
957 61756 : }
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 43 : DisplacedProblem::addCachedResidualDirectly(NumericVector<Number> & residual, const THREAD_ID tid)
968 : {
969 86 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
970 43 : _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 86 : if (_displaced_solver_systems[currentNlSysNum()]->hasVector(
977 43 : _displaced_solver_systems[currentNlSysNum()]->nonTimeVectorTag()))
978 129 : _assembly[tid][currentNlSysNum()]->addCachedResidualDirectly(
979 : residual,
980 86 : Assembly::GlobalDataKey{},
981 43 : 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 43 : _assembly[tid][currentNlSysNum()]->clearCachedResiduals(Assembly::GlobalDataKey{});
986 43 : }
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 192 : DisplacedProblem::addJacobian(const THREAD_ID tid)
1008 : {
1009 192 : _assembly[tid][currentNlSysNum()]->addJacobian(Assembly::GlobalDataKey{});
1010 192 : }
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 50 : DisplacedProblem::addJacobianNeighbor(const THREAD_ID tid)
1020 : {
1021 50 : _assembly[tid][currentNlSysNum()]->addJacobianNeighbor(Assembly::GlobalDataKey{});
1022 50 : }
1023 :
1024 : void
1025 3072 : DisplacedProblem::addJacobianNeighborLowerD(const THREAD_ID tid)
1026 : {
1027 3072 : _assembly[tid][currentNlSysNum()]->addJacobianNeighborLowerD(Assembly::GlobalDataKey{});
1028 3072 : }
1029 :
1030 : void
1031 192 : DisplacedProblem::addJacobianLowerD(const THREAD_ID tid)
1032 : {
1033 192 : _assembly[tid][currentNlSysNum()]->addJacobianLowerD(Assembly::GlobalDataKey{});
1034 192 : }
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 745052 : DisplacedProblem::prepareShapes(unsigned int var, const THREAD_ID tid)
1091 : {
1092 745052 : _assembly[tid][currentNlSysNum()]->copyShapes(var);
1093 745052 : }
1094 :
1095 : void
1096 8872 : DisplacedProblem::prepareFaceShapes(unsigned int var, const THREAD_ID tid)
1097 : {
1098 8872 : _assembly[tid][currentNlSysNum()]->copyFaceShapes(var);
1099 8872 : }
1100 :
1101 : void
1102 3352 : DisplacedProblem::prepareNeighborShapes(unsigned int var, const THREAD_ID tid)
1103 : {
1104 3352 : _assembly[tid][currentNlSysNum()]->copyNeighborShapes(var);
1105 3352 : }
1106 :
1107 : void
1108 4135 : DisplacedProblem::updateGeomSearch(GeometricSearchData::GeometricSearchType type)
1109 : {
1110 4135 : TIME_SECTION("updateGeometricSearch", 3, "Updating Displaced GeometricSearch");
1111 :
1112 4135 : _geometric_search_data.update(type);
1113 4135 : }
1114 :
1115 : void
1116 509 : 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 509 : _eq.reinit();
1122 509 : if (contract_mesh)
1123 : // Once vectors are restricted, we can delete children of coarsened elements
1124 454 : _mesh.getMesh().contract();
1125 509 : 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 454 : MeshRefinement refinement(_mesh.getMesh());
1130 454 : refinement.clean_refinement_flags();
1131 454 : }
1132 :
1133 : // Since the mesh has changed, we need to make sure that we update any of our
1134 : // MOOSE-system specific data.
1135 1018 : for (auto & nl : _displaced_solver_systems)
1136 509 : nl->reinit();
1137 509 : _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 509 : _mesh.clearQuadratureNodes();
1143 :
1144 509 : _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 509 : updateMesh(/*mesh_changing=*/true);
1150 509 : }
1151 :
1152 : void
1153 805257 : DisplacedProblem::addGhostedElem(dof_id_type elem_id)
1154 : {
1155 805257 : _mproblem.addGhostedElem(elem_id);
1156 805257 : }
1157 :
1158 : void
1159 26710 : DisplacedProblem::addGhostedBoundary(BoundaryID boundary_id)
1160 : {
1161 26710 : _mproblem.addGhostedBoundary(boundary_id);
1162 26710 : }
1163 :
1164 : void
1165 0 : DisplacedProblem::ghostGhostedBoundaries()
1166 : {
1167 0 : _mproblem.ghostGhostedBoundaries();
1168 0 : }
1169 :
1170 : MooseMesh &
1171 416084 : DisplacedProblem::refMesh()
1172 : {
1173 416084 : 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 579 : 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 1158 : NodeRange node_range(_mesh.getMesh().nodes_begin(),
1219 1158 : _mesh.getMesh().nodes_end(),
1220 579 : /*grainsize=*/1);
1221 :
1222 579 : ResetDisplacedMeshThread rdmt(_mproblem, *this);
1223 :
1224 : // Undisplace the mesh using threads.
1225 579 : Threads::parallel_reduce(node_range, rdmt);
1226 579 : }
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 501024 : DisplacedProblem::computingScalingJacobian() const
1242 : {
1243 501024 : return _mproblem.computingScalingJacobian();
1244 : }
1245 :
1246 : bool
1247 0 : DisplacedProblem::computingScalingResidual() const
1248 : {
1249 0 : return _mproblem.computingScalingResidual();
1250 : }
1251 :
1252 : void
1253 2051 : DisplacedProblem::initialSetup()
1254 : {
1255 2051 : SubProblem::initialSetup();
1256 :
1257 4102 : for (auto & nl : _displaced_solver_systems)
1258 2051 : nl->initialSetup();
1259 2051 : _displaced_aux->initialSetup();
1260 2051 : }
1261 :
1262 : void
1263 30557 : DisplacedProblem::timestepSetup()
1264 : {
1265 30557 : SubProblem::timestepSetup();
1266 :
1267 61114 : for (auto & nl : _displaced_solver_systems)
1268 30557 : nl->timestepSetup();
1269 30557 : _displaced_aux->timestepSetup();
1270 30557 : }
1271 :
1272 : void
1273 144802 : DisplacedProblem::customSetup(const ExecFlagType & exec_type)
1274 : {
1275 144802 : SubProblem::customSetup(exec_type);
1276 :
1277 289604 : for (auto & nl : _displaced_solver_systems)
1278 144802 : nl->customSetup(exec_type);
1279 144802 : _displaced_aux->customSetup(exec_type);
1280 144802 : }
1281 :
1282 : void
1283 126628 : DisplacedProblem::residualSetup()
1284 : {
1285 126628 : SubProblem::residualSetup();
1286 :
1287 253256 : for (auto & nl : _displaced_solver_systems)
1288 126628 : nl->residualSetup();
1289 126628 : _displaced_aux->residualSetup();
1290 126628 : }
1291 :
1292 : void
1293 21200 : DisplacedProblem::jacobianSetup()
1294 : {
1295 21200 : SubProblem::jacobianSetup();
1296 :
1297 42400 : for (auto & nl : _displaced_solver_systems)
1298 21200 : nl->jacobianSetup();
1299 21200 : _displaced_aux->jacobianSetup();
1300 21200 : }
1301 :
1302 : void
1303 288 : DisplacedProblem::haveADObjects(const bool have_ad_objects)
1304 : {
1305 288 : _have_ad_objects = have_ad_objects;
1306 288 : _mproblem.SubProblem::haveADObjects(have_ad_objects);
1307 288 : }
1308 :
1309 : std::pair<bool, unsigned int>
1310 32547 : DisplacedProblem::determineSolverSystem(const std::string & var_name,
1311 : const bool error_if_not_found) const
1312 : {
1313 32547 : return _mproblem.determineSolverSystem(var_name, error_if_not_found);
1314 : }
1315 :
1316 : Assembly &
1317 53234214 : 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 53234214 : return *_assembly[tid][sys_num];
1323 : }
1324 :
1325 : const Assembly &
1326 45133 : 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 45133 : return *_assembly[tid][sys_num];
1332 : }
1333 :
1334 : std::size_t
1335 9351090 : DisplacedProblem::numNonlinearSystems() const
1336 : {
1337 9351090 : 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 362956 : DisplacedProblem::numSolverSystems() const
1348 : {
1349 362956 : return _mproblem.numSolverSystems();
1350 : }
1351 :
1352 : const std::vector<VectorTag> &
1353 8306035 : DisplacedProblem::currentResidualVectorTags() const
1354 : {
1355 8306035 : return _mproblem.currentResidualVectorTags();
1356 : }
1357 :
1358 : bool
1359 54856892 : DisplacedProblem::safeAccessTaggedMatrices() const
1360 : {
1361 54856892 : return _mproblem.safeAccessTaggedMatrices();
1362 : }
1363 :
1364 : bool
1365 441 : DisplacedProblem::safeAccessTaggedVectors() const
1366 : {
1367 441 : return _mproblem.safeAccessTaggedVectors();
1368 : }
1369 :
1370 : void
1371 0 : DisplacedProblem::needFV()
1372 : {
1373 0 : _mproblem.needFV();
1374 0 : }
1375 :
1376 : bool
1377 300471 : DisplacedProblem::haveFV() const
1378 : {
1379 300471 : return _mproblem.haveFV();
1380 : }
1381 :
1382 : bool
1383 1592366 : DisplacedProblem::hasNonlocalCoupling() const
1384 : {
1385 1592366 : 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 2262 : DisplacedProblem::nonlocalCouplingMatrix(const unsigned i) const
1408 : {
1409 2262 : return _mproblem.nonlocalCouplingMatrix(i);
1410 : }
1411 :
1412 : bool
1413 0 : DisplacedProblem::checkNonlocalCouplingRequirement() const
1414 : {
1415 0 : return _mproblem.checkNonlocalCouplingRequirement();
1416 : }
1417 :
1418 150580 : DisplacedProblem::UpdateDisplacedMeshThread::UpdateDisplacedMeshThread(
1419 150580 : FEProblemBase & fe_problem, DisplacedProblem & displaced_problem)
1420 : : ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>(fe_problem),
1421 150580 : _displaced_problem(displaced_problem),
1422 301160 : _ref_mesh(_displaced_problem.refMesh()),
1423 150580 : _nl_soln(_displaced_problem._nl_solution),
1424 150580 : _aux_soln(*_displaced_problem._aux_solution),
1425 150580 : _has_displacement(false)
1426 : {
1427 150580 : this->init();
1428 150580 : }
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 150580 : DisplacedProblem::UpdateDisplacedMeshThread::init()
1445 : {
1446 150580 : std::vector<std::string> & displacement_variables = _displaced_problem._displacements;
1447 150580 : unsigned int num_displacements = displacement_variables.size();
1448 150580 : auto & es = _displaced_problem.es();
1449 :
1450 150580 : _sys_to_var_num_and_direction.clear();
1451 150580 : _sys_to_nonghost_and_ghost_soln.clear();
1452 :
1453 502360 : for (unsigned int i = 0; i < num_displacements; i++)
1454 : {
1455 351780 : std::string displacement_name = displacement_variables[i];
1456 :
1457 459855 : for (const auto sys_num : make_range(es.n_systems()))
1458 : {
1459 459855 : auto & sys = es.get_system(sys_num);
1460 459855 : if (sys.has_variable(displacement_name))
1461 : {
1462 351780 : auto & val = _sys_to_var_num_and_direction[sys.number()];
1463 351780 : val.first.push_back(sys.variable_number(displacement_name));
1464 351780 : val.second.push_back(i);
1465 351780 : break;
1466 : }
1467 : }
1468 351780 : }
1469 :
1470 303364 : for (const auto & pr : _sys_to_var_num_and_direction)
1471 : {
1472 152784 : 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 152784 : sys.number() < _nl_soln.size() ? _nl_soln[sys.number()] : &_aux_soln;
1478 305568 : _sys_to_nonghost_and_ghost_soln.emplace(
1479 152784 : sys.number(),
1480 0 : std::make_pair(nonghost_soln,
1481 305568 : NumericVector<Number>::build(nonghost_soln->comm()).release()));
1482 : }
1483 :
1484 150580 : ConstNodeRange node_range(_ref_mesh.getMesh().nodes_begin(), _ref_mesh.getMesh().nodes_end());
1485 :
1486 303364 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1487 : {
1488 152784 : auto & sys = es.get_system(sys_num);
1489 : AllNodesSendListThread send_list(
1490 152784 : this->_fe_problem, _ref_mesh, var_num_and_direction.first, sys);
1491 152784 : Threads::parallel_reduce(node_range, send_list);
1492 152784 : send_list.unique();
1493 152784 : auto & [soln, ghost_soln] = libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num);
1494 305568 : ghost_soln->init(
1495 152784 : soln->size(), soln->local_size(), send_list.send_list(), true, libMesh::GHOSTED);
1496 152784 : soln->localize(*ghost_soln, send_list.send_list());
1497 152784 : }
1498 :
1499 150580 : _has_displacement = false;
1500 150580 : }
1501 :
1502 : void
1503 18225870 : DisplacedProblem::UpdateDisplacedMeshThread::onNode(NodeRange::const_iterator & nd)
1504 : {
1505 18225870 : Node & displaced_node = *(*nd);
1506 :
1507 18225870 : Node & reference_node = _ref_mesh.nodeRef(displaced_node.id());
1508 :
1509 36902740 : for (auto & [sys_num, var_num_and_direction] : _sys_to_var_num_and_direction)
1510 : {
1511 18676870 : auto & var_numbers = var_num_and_direction.first;
1512 18676870 : auto & directions = var_num_and_direction.second;
1513 57436126 : for (const auto i : index_range(var_numbers))
1514 : {
1515 38759256 : const auto direction = directions[i];
1516 38759256 : if (reference_node.n_dofs(sys_num, var_numbers[i]) > 0)
1517 : {
1518 38759128 : Real coord = reference_node(direction) +
1519 77518256 : (*libmesh_map_find(_sys_to_nonghost_and_ghost_soln, sys_num).second)(
1520 38759128 : reference_node.dof_number(sys_num, var_numbers[i], 0));
1521 38759128 : if (displaced_node(direction) != coord)
1522 : {
1523 3267632 : displaced_node(direction) = coord;
1524 3267632 : _has_displacement = true;
1525 : }
1526 : }
1527 : }
1528 : }
1529 18225870 : }
|
