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 : #pragma once
11 :
12 : #include "SubProblem.h"
13 : #include "DisplacedSystem.h"
14 : #include "GeometricSearchData.h"
15 : #include "ThreadedNodeLoop.h"
16 :
17 : // libMesh
18 : #include "libmesh/equation_systems.h"
19 : #include "libmesh/enum_quadrature_type.h"
20 :
21 : // Forward declarations
22 : class MooseVariableFieldBase;
23 : class AssemblyData;
24 : class MooseMesh;
25 : class Assembly;
26 : class FEProblemBase;
27 : class LineSearch;
28 :
29 : // libMesh forward declarations
30 : namespace libMesh
31 : {
32 : template <typename T>
33 : class NumericVector;
34 : class DofMap;
35 : }
36 :
37 : class DisplacedProblem : public SubProblem
38 : {
39 : public:
40 : static InputParameters validParams();
41 :
42 : DisplacedProblem(DisplacedProblem &&) = delete;
43 : DisplacedProblem & operator=(DisplacedProblem &&) = delete;
44 :
45 : DisplacedProblem(const InputParameters & parameters);
46 : ~DisplacedProblem();
47 :
48 165737 : virtual EquationSystems & es() override { return _eq; }
49 414085 : virtual MooseMesh & mesh() override { return _mesh; }
50 43714386 : virtual const MooseMesh & mesh() const override { return _mesh; }
51 0 : const MooseMesh & mesh(bool libmesh_dbg_var(use_displaced)) const override
52 : {
53 : mooseAssert(use_displaced, "An undisplaced mesh was queried from the displaced problem");
54 0 : return mesh();
55 : }
56 : MooseMesh & refMesh();
57 :
58 : DisplacedSystem & solverSys(const unsigned int sys_num);
59 48664 : DisplacedSystem & auxSys() { return *_displaced_aux; }
60 :
61 : virtual const SystemBase & systemBaseNonlinear(const unsigned int sys_num) const override;
62 : virtual SystemBase & systemBaseNonlinear(const unsigned int sys_num) override;
63 :
64 : virtual const SystemBase & systemBaseLinear(const unsigned int sys_num) const override;
65 : virtual SystemBase & systemBaseLinear(const unsigned int sys_num) override;
66 :
67 : virtual const SystemBase & systemBaseSolver(const unsigned int sys_num) const override;
68 : virtual SystemBase & systemBaseSolver(const unsigned int sys_num) override;
69 :
70 0 : virtual const SystemBase & systemBaseAuxiliary() const override { return *_displaced_aux; }
71 522163 : virtual SystemBase & systemBaseAuxiliary() override { return *_displaced_aux; }
72 :
73 : // Return a constant reference to the vector of variable names.
74 339 : const std::vector<std::string> & getDisplacementVarNames() const { return _displacements; }
75 :
76 : virtual void createQRules(QuadratureType type,
77 : Order order,
78 : Order volume_order,
79 : Order face_order,
80 : SubdomainID block,
81 : bool allow_negative_qweights = true);
82 :
83 : void bumpVolumeQRuleOrder(Order order, SubdomainID block);
84 : void bumpAllQRuleOrder(Order order, SubdomainID block);
85 :
86 : virtual void init() override;
87 : virtual bool solverSystemConverged(const unsigned int solver_sys_num) override;
88 : virtual unsigned int nlSysNum(const NonlinearSystemName & nl_sys_name) const override;
89 : virtual unsigned int linearSysNum(const LinearSystemName & sys_name) const override;
90 : virtual unsigned int solverSysNum(const SolverSystemName & sys_name) const override;
91 :
92 : /// Get the time integrators from the problem
93 : void addTimeIntegrator();
94 :
95 : /**
96 : * Allocate vectors and save old solutions into them.
97 : */
98 : virtual void saveOldSolutions();
99 :
100 : /**
101 : * Restore old solutions from the backup vectors and deallocate them.
102 : */
103 : virtual void restoreOldSolutions();
104 :
105 : /**
106 : * Copy the provided solution into the displaced auxiliary system
107 : */
108 : void syncAuxSolution(const NumericVector<Number> & aux_soln);
109 :
110 : /**
111 : * Copy the solutions on the undisplaced systems to the displaced systems.
112 : */
113 : void syncSolutions();
114 :
115 : /**
116 : * Synchronize the solutions on the displaced systems to the given solutions. The nonlinear
117 : * solutions argument is a map from the nonlinear system number to the solution that we want to
118 : * set that nonlinear system's solution to
119 : */
120 : void syncSolutions(const std::map<unsigned int, const NumericVector<Number> *> & nl_solns,
121 : const NumericVector<Number> & aux_soln);
122 :
123 : /**
124 : * Copy the solutions on the undisplaced systems to the displaced systems and
125 : * reinitialize the geometry search data and Dirac kernel information due to mesh displacement.
126 : * The parameter \p mesh_changing indicates whether this method is getting called because of mesh
127 : * changes, e.g. due to mesh adaptivity. If \p mesh_changing we need to renitialize the
128 : * GeometricSearchData instead of simply update. Reinitialization operations are a super-set of
129 : * update operations. Reinitialization for example re-generates neighbor nodes in
130 : * NearestNodeLocators, while update does not. Additionally we do not want to use the undisplaced
131 : * mesh solution because it may be out-of-sync, whereas our displaced mesh solution should be in
132 : * the correct state after getting restricted/prolonged in EquationSystems::reinit
133 : */
134 : virtual void updateMesh(bool mesh_changing = false);
135 :
136 : /**
137 : * Synchronize the solutions on the displaced systems to the given solutions and
138 : * reinitialize the geometry search data and Dirac kernel information due to mesh displacement.
139 : */
140 : virtual void updateMesh(const std::map<unsigned int, const NumericVector<Number> *> & nl_soln,
141 : const NumericVector<Number> & aux_soln);
142 :
143 : virtual TagID addVectorTag(const TagName & tag_name,
144 : const Moose::VectorTagType type = Moose::VECTOR_TAG_RESIDUAL) override;
145 : virtual const VectorTag & getVectorTag(const TagID tag_id) const override;
146 : virtual TagID getVectorTagID(const TagName & tag_name) const override;
147 : virtual TagName vectorTagName(const TagID tag_id) const override;
148 : virtual bool vectorTagExists(const TagID tag_id) const override;
149 : virtual bool vectorTagExists(const TagName & tag_name) const override;
150 : virtual unsigned int
151 : numVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const override;
152 : virtual const std::vector<VectorTag> &
153 : getVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const override;
154 : virtual Moose::VectorTagType vectorTagType(const TagID tag_id) const override;
155 :
156 : virtual TagID addMatrixTag(TagName tag_name) override;
157 : virtual TagID getMatrixTagID(const TagName & tag_name) const override;
158 : virtual TagName matrixTagName(TagID tag) override;
159 : virtual bool matrixTagExists(const TagName & tag_name) const override;
160 : virtual bool matrixTagExists(TagID tag_id) const override;
161 : virtual unsigned int numMatrixTags() const override;
162 : virtual bool safeAccessTaggedMatrices() const override;
163 : virtual bool safeAccessTaggedVectors() const override;
164 :
165 : virtual bool isTransient() const override;
166 :
167 : // Variables /////
168 : virtual bool hasVariable(const std::string & var_name) const override;
169 : using SubProblem::getVariable;
170 : virtual const MooseVariableFieldBase &
171 : getVariable(const THREAD_ID tid,
172 : const std::string & var_name,
173 : Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
174 : Moose::VarFieldType expected_var_field_type =
175 : Moose::VarFieldType::VAR_FIELD_ANY) const override;
176 : virtual MooseVariable & getStandardVariable(const THREAD_ID tid,
177 : const std::string & var_name) override;
178 : virtual MooseVariableFieldBase & getActualFieldVariable(const THREAD_ID tid,
179 : const std::string & var_name) override;
180 : virtual VectorMooseVariable & getVectorVariable(const THREAD_ID tid,
181 : const std::string & var_name) override;
182 : virtual ArrayMooseVariable & getArrayVariable(const THREAD_ID tid,
183 : const std::string & var_name) override;
184 : virtual bool hasScalarVariable(const std::string & var_name) const override;
185 : virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid,
186 : const std::string & var_name) override;
187 : virtual System & getSystem(const std::string & var_name) override;
188 :
189 : virtual void addVariable(const std::string & var_type,
190 : const std::string & name,
191 : InputParameters & parameters,
192 : unsigned int nl_system_number);
193 : virtual void addAuxVariable(const std::string & var_type,
194 : const std::string & name,
195 : InputParameters & parameters);
196 : //
197 : // Adaptivity /////
198 : virtual void initAdaptivity();
199 : void meshChanged(bool contract_mesh, bool clean_refinement_flags);
200 :
201 : // reinit /////
202 : virtual void prepare(const Elem * elem, const THREAD_ID tid) override;
203 : virtual void prepareNonlocal(const THREAD_ID tid);
204 : virtual void prepareFace(const Elem * elem, const THREAD_ID tid) override;
205 : virtual void prepare(const Elem * elem,
206 : unsigned int ivar,
207 : unsigned int jvar,
208 : const std::vector<dof_id_type> & dof_indices,
209 : const THREAD_ID tid) override;
210 : virtual void setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid) override;
211 : virtual void
212 : setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
213 : virtual void prepareBlockNonlocal(unsigned int ivar,
214 : unsigned int jvar,
215 : const std::vector<dof_id_type> & idof_indices,
216 : const std::vector<dof_id_type> & jdof_indices,
217 : const THREAD_ID tid);
218 : virtual void prepareAssembly(const THREAD_ID tid) override;
219 : virtual void prepareAssemblyNeighbor(const THREAD_ID tid);
220 :
221 : virtual bool reinitDirac(const Elem * elem, const THREAD_ID tid) override;
222 :
223 : virtual void reinitElem(const Elem * elem, const THREAD_ID tid) override;
224 : virtual void reinitElemPhys(const Elem * elem,
225 : const std::vector<Point> & phys_points_in_elem,
226 : const THREAD_ID tid) override;
227 : virtual void reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
228 : virtual void reinitNode(const Node * node, const THREAD_ID tid) override;
229 : virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid) override;
230 : virtual void reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) override;
231 : virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes,
232 : const THREAD_ID tid) override;
233 : virtual void reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
234 :
235 : /**
236 : * reinitialize neighbor routine
237 : * @param elem The element driving the reinit (note that this is not the *neighbor*)
238 : * @param side The side, e.g. face, of the \p elem that we want to reinit
239 : * @param tid The thread for which we are reiniting
240 : * @param neighbor_reference_points Specify the referrence points for the
241 : * neighbor element. Useful if the element and neighbor faces are
242 : * not coincident
243 : */
244 : void reinitNeighbor(const Elem * elem,
245 : unsigned int side,
246 : const THREAD_ID tid,
247 : const std::vector<Point> * neighbor_reference_points);
248 :
249 : virtual void reinitNeighborPhys(const Elem * neighbor,
250 : unsigned int neighbor_side,
251 : const std::vector<Point> & physical_points,
252 : const THREAD_ID tid) override;
253 : virtual void reinitNeighborPhys(const Elem * neighbor,
254 : const std::vector<Point> & physical_points,
255 : const THREAD_ID tid) override;
256 : virtual void
257 : reinitElemNeighborAndLowerD(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
258 : virtual void reinitScalars(const THREAD_ID tid,
259 : bool reinit_for_derivative_reordering = false) override;
260 : virtual void reinitOffDiagScalars(const THREAD_ID tid) override;
261 :
262 : /// Fills "elems" with the elements that should be looped over for Dirac Kernels
263 : virtual void getDiracElements(std::set<const Elem *> & elems) override;
264 : virtual void clearDiracInfo() override;
265 :
266 : virtual void addResidual(const THREAD_ID tid) override;
267 : virtual void addResidualNeighbor(const THREAD_ID tid) override;
268 : virtual void addResidualLower(const THREAD_ID tid) override;
269 :
270 : virtual void addCachedResidualDirectly(NumericVector<Number> & residual, const THREAD_ID tid);
271 :
272 : virtual void setResidual(NumericVector<Number> & residual, const THREAD_ID tid) override;
273 : virtual void setResidualNeighbor(NumericVector<Number> & residual, const THREAD_ID tid) override;
274 :
275 : virtual void addJacobian(const THREAD_ID tid) override;
276 : virtual void addJacobianNonlocal(const THREAD_ID tid);
277 : virtual void addJacobianNeighbor(const THREAD_ID tid) override;
278 : virtual void addJacobianNeighborLowerD(const THREAD_ID tid) override;
279 : virtual void addJacobianLowerD(const THREAD_ID tid) override;
280 : virtual void addJacobianBlockTags(SparseMatrix<Number> & jacobian,
281 : unsigned int ivar,
282 : unsigned int jvar,
283 : const libMesh::DofMap & dof_map,
284 : std::vector<dof_id_type> & dof_indices,
285 : const std::set<TagID> & tags,
286 : const THREAD_ID tid);
287 : void addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
288 : unsigned int ivar,
289 : unsigned int jvar,
290 : const libMesh::DofMap & dof_map,
291 : const std::vector<dof_id_type> & idof_indices,
292 : const std::vector<dof_id_type> & jdof_indices,
293 : const std::set<TagID> & tags,
294 : const THREAD_ID tid);
295 : virtual void addJacobianNeighbor(SparseMatrix<Number> & jacobian,
296 : unsigned int ivar,
297 : unsigned int jvar,
298 : const libMesh::DofMap & dof_map,
299 : std::vector<dof_id_type> & dof_indices,
300 : std::vector<dof_id_type> & neighbor_dof_indices,
301 : const std::set<TagID> & tags,
302 : const THREAD_ID tid) override;
303 :
304 : virtual void cacheJacobianNonlocal(const THREAD_ID tid);
305 :
306 : virtual void prepareShapes(unsigned int var, const THREAD_ID tid) override;
307 : virtual void prepareFaceShapes(unsigned int var, const THREAD_ID tid) override;
308 : virtual void prepareNeighborShapes(unsigned int var, const THREAD_ID tid) override;
309 :
310 : virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override;
311 : virtual const Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) const override;
312 :
313 : // Geom Search /////
314 : virtual void updateGeomSearch(
315 : GeometricSearchData::GeometricSearchType type = GeometricSearchData::ALL) override;
316 64311 : virtual GeometricSearchData & geomSearchData() override { return _geometric_search_data; }
317 :
318 : virtual bool computingPreSMOResidual(const unsigned int nl_sys_num) const override;
319 :
320 : virtual void onTimestepBegin() override;
321 : virtual void onTimestepEnd() override;
322 :
323 : /**
324 : * Return the list of elements that should have their DoFs ghosted to this processor.
325 : * @return The list
326 : */
327 : virtual std::set<dof_id_type> & ghostedElems() override;
328 :
329 : /**
330 : * Will make sure that all dofs connected to elem_id are ghosted to this processor
331 : */
332 : virtual void addGhostedElem(dof_id_type elem_id) override;
333 :
334 : /**
335 : * Will make sure that all necessary elements from boundary_id are ghosted to this processor
336 : * @param boundary_id Boundary ID
337 : */
338 : virtual void addGhostedBoundary(BoundaryID boundary_id) override;
339 :
340 : /**
341 : * Causes the boundaries added using addGhostedBoundary to actually be ghosted.
342 : */
343 : virtual void ghostGhostedBoundaries() override;
344 :
345 : /**
346 : * Resets the displaced mesh to the reference mesh. Required when
347 : * refining the DisplacedMesh.
348 : */
349 : void undisplaceMesh();
350 :
351 : virtual LineSearch * getLineSearch() override;
352 :
353 : virtual const CouplingMatrix * couplingMatrix(const unsigned int nl_sys_num) const override;
354 :
355 0 : virtual bool haveDisplaced() const override final { return true; }
356 :
357 : virtual bool computingScalingJacobian() const override final;
358 :
359 : virtual bool computingScalingResidual() const override final;
360 :
361 : virtual void initialSetup() override;
362 : virtual void timestepSetup() override;
363 : virtual void customSetup(const ExecFlagType & exec_type) override;
364 : virtual void residualSetup() override;
365 : virtual void jacobianSetup() override;
366 :
367 : using SubProblem::haveADObjects;
368 : virtual void haveADObjects(bool have_ad_objects) override;
369 :
370 : virtual std::size_t numNonlinearSystems() const override;
371 :
372 : virtual std::size_t numLinearSystems() const override;
373 :
374 : virtual std::size_t numSolverSystems() const override;
375 :
376 : virtual unsigned int currentNlSysNum() const override;
377 :
378 : virtual unsigned int currentLinearSysNum() const override;
379 :
380 : virtual const std::vector<VectorTag> & currentResidualVectorTags() const override;
381 :
382 : virtual void needFV() override;
383 : virtual bool haveFV() const override;
384 :
385 : virtual bool hasNonlocalCoupling() const override;
386 : virtual bool checkNonlocalCouplingRequirement() const override;
387 : virtual const libMesh::CouplingMatrix & nonlocalCouplingMatrix(const unsigned i) const override;
388 :
389 : protected:
390 : FEProblemBase & _mproblem;
391 : MooseMesh & _mesh;
392 : EquationSystems _eq;
393 : /// reference mesh
394 : MooseMesh & _ref_mesh;
395 : std::vector<std::string> _displacements;
396 :
397 : std::vector<std::unique_ptr<DisplacedSystem>> _displaced_solver_systems;
398 : std::unique_ptr<DisplacedSystem> _displaced_aux;
399 :
400 : /// The nonlinear system solutions
401 : std::vector<const NumericVector<Number> *> _nl_solution;
402 :
403 : /// The auxiliary system solution
404 : const NumericVector<Number> * _aux_solution;
405 :
406 : std::vector<std::vector<std::unique_ptr<Assembly>>> _assembly;
407 :
408 : GeometricSearchData _geometric_search_data;
409 :
410 : class UpdateDisplacedMeshThread : public ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>
411 : {
412 : public:
413 : UpdateDisplacedMeshThread(FEProblemBase & fe_problem, DisplacedProblem & displaced_problem);
414 :
415 : UpdateDisplacedMeshThread(UpdateDisplacedMeshThread & x, Threads::split split);
416 :
417 : virtual void onNode(NodeRange::const_iterator & nd) override;
418 :
419 14492 : void join(const UpdateDisplacedMeshThread & y)
420 : {
421 14492 : if (y._has_displacement)
422 5387 : _has_displacement = true;
423 14492 : }
424 :
425 : /**
426 : * Whether the displaced mesh is modified by the latest call to operator()
427 : */
428 150576 : bool hasDisplacement()
429 : {
430 : mooseAssert(!Threads::in_threads,
431 : "This function requires a MPI all-gathering operation that cannot be in a "
432 : "threaded scope.");
433 150576 : _ref_mesh.comm().max(_has_displacement);
434 150576 : return _has_displacement;
435 : }
436 :
437 : protected:
438 : void init();
439 :
440 : /// Diplaced problem
441 : DisplacedProblem & _displaced_problem;
442 : /// Original mesh
443 : MooseMesh & _ref_mesh;
444 : /// Solution vectors of the nonlinear systems on the displaced problem
445 : const std::vector<const NumericVector<Number> *> & _nl_soln;
446 : /// Solution vector of the auxliary system on the displaced problem
447 : const NumericVector<Number> & _aux_soln;
448 :
449 : // Solution vectors with expanded ghosting, for ReplicatedMesh or
450 : // for DistributedMesh cases where the standard algebraic ghosting
451 : // doesn't reach as far as the geometric ghosting
452 : std::map<unsigned int,
453 : std::pair<const NumericVector<Number> *, std::shared_ptr<NumericVector<Number>>>>
454 : _sys_to_nonghost_and_ghost_soln;
455 :
456 : private:
457 : /// To locate the system numbers, variable numbers of all displacement variables
458 : std::map<unsigned int, std::pair<std::vector<unsigned int>, std::vector<unsigned int>>>
459 : _sys_to_var_num_and_direction;
460 :
461 : /// A flag to be set by operator() for indicating whether the displaced mesh is
462 : /// indeed modified
463 : bool _has_displacement;
464 : };
465 :
466 : private:
467 : virtual std::pair<bool, unsigned int>
468 : determineSolverSystem(const std::string & var_name,
469 : bool error_if_not_found = false) const override;
470 :
471 : friend class UpdateDisplacedMeshThread;
472 : friend class Restartable;
473 : };
474 :
475 : inline DisplacedSystem &
476 59892 : DisplacedProblem::solverSys(const unsigned int sys_num)
477 : {
478 : mooseAssert(sys_num < _displaced_solver_systems.size(),
479 : "System number greater than the number of nonlinear systems");
480 59892 : return *_displaced_solver_systems[sys_num];
481 : }
482 :
483 : inline const SystemBase &
484 516 : DisplacedProblem::systemBaseNonlinear(const unsigned int sys_num) const
485 : {
486 : mooseAssert(sys_num < _displaced_solver_systems.size(),
487 : "System number greater than the number of nonlinear systems");
488 516 : return *_displaced_solver_systems[sys_num];
489 : }
490 :
491 : inline SystemBase &
492 99967 : DisplacedProblem::systemBaseNonlinear(const unsigned int sys_num)
493 : {
494 : mooseAssert(sys_num < _displaced_solver_systems.size(),
495 : "System number greater than the number of nonlinear systems");
496 99967 : return *_displaced_solver_systems[sys_num];
497 : }
498 :
499 : inline const SystemBase &
500 0 : DisplacedProblem::systemBaseLinear(const unsigned int /*sys_num*/) const
501 : {
502 0 : mooseError("Linear systems are not supported for displaced problems yet.");
503 : }
504 :
505 : inline SystemBase &
506 0 : DisplacedProblem::systemBaseLinear(const unsigned int /*sys_num*/)
507 : {
508 0 : mooseError("Linear systems are not supported for displaced problems yet.");
509 : }
510 :
511 : inline const SystemBase &
512 0 : DisplacedProblem::systemBaseSolver(const unsigned int sys_num) const
513 : {
514 : mooseAssert(sys_num < _displaced_solver_systems.size(),
515 : "System number greater than the number of solver systems");
516 0 : return *_displaced_solver_systems[sys_num];
517 : }
518 :
519 : inline SystemBase &
520 362956 : DisplacedProblem::systemBaseSolver(const unsigned int sys_num)
521 : {
522 : mooseAssert(sys_num < _displaced_solver_systems.size(),
523 : "System number greater than the number of solver systems");
524 362956 : return *_displaced_solver_systems[sys_num];
525 : }
|
