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 : #include "MooseVariableScalar.h"
11 : #include "SubProblem.h"
12 : #include "SystemBase.h"
13 : #include "Assembly.h"
14 : #include "SystemBase.h"
15 :
16 : // libMesh
17 : #include "libmesh/numeric_vector.h"
18 : #include "libmesh/dof_map.h"
19 :
20 : #include <limits>
21 :
22 : registerMooseObject("MooseApp", MooseVariableScalar);
23 :
24 : InputParameters
25 6353 : MooseVariableScalar::validParams()
26 : {
27 6353 : auto params = MooseVariableBase::validParams();
28 12706 : params.addClassDescription("Moose wrapper class around scalar variables");
29 19059 : params.set<MooseEnum>("family") = "SCALAR";
30 6353 : return params;
31 0 : }
32 :
33 3136 : MooseVariableScalar::MooseVariableScalar(const InputParameters & parameters)
34 : : MooseVariableBase(parameters),
35 : Moose::FunctorBase<ADReal>(name()),
36 3136 : _need_u_dot(false),
37 3136 : _need_u_dotdot(false),
38 3136 : _need_u_dot_old(false),
39 3136 : _need_u_dotdot_old(false),
40 3136 : _need_du_dot_du(false),
41 3136 : _need_du_dotdot_du(false),
42 3136 : _need_u_old(false),
43 3136 : _need_u_older(false),
44 3136 : _need_ad(false),
45 3136 : _need_ad_u(false),
46 12544 : _need_ad_u_dot(false)
47 : {
48 3136 : auto num_vector_tags = _sys.subproblem().numVectorTags();
49 :
50 3136 : _vector_tag_u.resize(num_vector_tags);
51 3136 : _need_vector_tag_u.resize(num_vector_tags);
52 6272 : }
53 :
54 : void
55 3124 : MooseVariableScalar::sizeMatrixTagData()
56 : {
57 3124 : auto num_matrix_tags = _sys.subproblem().numMatrixTags();
58 :
59 3124 : _matrix_tag_u.resize(num_matrix_tags);
60 3124 : _need_matrix_tag_u.resize(num_matrix_tags);
61 3124 : }
62 :
63 2984 : MooseVariableScalar::~MooseVariableScalar()
64 : {
65 2984 : _u.release();
66 2984 : _u_old.release();
67 2984 : _u_older.release();
68 :
69 2984 : _u_dot.release();
70 2984 : _u_dotdot.release();
71 2984 : _u_dot_old.release();
72 2984 : _u_dotdot_old.release();
73 2984 : _du_dot_du.release();
74 2984 : _du_dotdot_du.release();
75 :
76 17344 : for (auto & _tag_u : _vector_tag_u)
77 14360 : _tag_u.release();
78 :
79 2984 : _vector_tag_u.clear();
80 :
81 9126 : for (auto & _tag_u : _matrix_tag_u)
82 6142 : _tag_u.release();
83 :
84 2984 : _matrix_tag_u.clear();
85 2984 : }
86 :
87 : void
88 482384 : MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/)
89 : {
90 482384 : if (reinit_for_derivative_reordering)
91 : {
92 : // We've already calculated everything in an earlier reinit. All we have to do is re-sort the
93 : // derivative vector for AD calculations
94 17860 : if (_need_ad)
95 593 : computeAD(/*nodal_ordering=*/true);
96 17860 : return;
97 : }
98 :
99 : // We want to make sure that we're not allocating data with any of the
100 : // accessors that follow, but _sys is non-const
101 464524 : const SystemBase & sys = _sys;
102 :
103 464524 : const NumericVector<Real> & current_solution = *sys.currentSolution();
104 464524 : const NumericVector<Real> * u_dot = sys.solutionUDot();
105 464524 : const NumericVector<Real> * u_dotdot = sys.solutionUDotDot();
106 464524 : const NumericVector<Real> * u_dot_old = sys.solutionUDotOld();
107 464524 : const NumericVector<Real> * u_dotdot_old = sys.solutionUDotDotOld();
108 464524 : const Real & du_dot_du = sys.duDotDu(number());
109 464524 : const Real & du_dotdot_du = sys.duDotDotDu();
110 464524 : auto safe_access_tagged_vectors = sys.subproblem().safeAccessTaggedVectors();
111 464524 : auto safe_access_tagged_matrices = sys.subproblem().safeAccessTaggedMatrices();
112 : auto & active_coupleable_matrix_tags =
113 464524 : sys.subproblem().getActiveScalarVariableCoupleableMatrixTags(_tid);
114 :
115 464524 : _dof_map.SCALAR_dof_indices(_dof_indices, _var_num);
116 :
117 464524 : auto n = _dof_indices.size();
118 464524 : _u.resize(n);
119 464524 : if (_need_u_old)
120 8186 : _u_old.resize(n);
121 464524 : if (_need_u_older)
122 2148 : _u_older.resize(n);
123 :
124 3562139 : for (auto & _tag_u : _vector_tag_u)
125 3097615 : _tag_u.resize(n);
126 :
127 1404418 : for (auto & _tag_u : _matrix_tag_u)
128 939894 : _tag_u.resize(n);
129 :
130 464524 : _du_dot_du.clear();
131 464524 : _du_dot_du.resize(n, du_dot_du);
132 :
133 464524 : if (_need_u_dot)
134 254012 : _u_dot.resize(n);
135 :
136 464524 : if (_need_u_dotdot)
137 0 : _u_dotdot.resize(n);
138 :
139 464524 : if (_need_u_dot_old)
140 0 : _u_dot_old.resize(n);
141 :
142 464524 : if (_need_u_dotdot_old)
143 0 : _u_dotdot_old.resize(n);
144 :
145 464524 : if (_need_du_dot_du)
146 : {
147 252223 : _du_dot_du.clear();
148 252223 : _du_dot_du.resize(n, du_dot_du);
149 : }
150 :
151 464524 : if (_need_du_dotdot_du)
152 : {
153 0 : _du_dotdot_du.clear();
154 0 : _du_dotdot_du.resize(n, du_dotdot_du);
155 : }
156 :
157 : // If we have an empty partition, or if we have a partition which
158 : // does not include any of the subdomains of a subdomain-restricted
159 : // variable, then we do not have access to that variable! Hopefully
160 : // we won't need the indices we lack.
161 464524 : if (_dof_map.all_semilocal_indices(_dof_indices))
162 : {
163 464524 : current_solution.get(_dof_indices, &_u[0]);
164 464524 : if (_need_u_old)
165 8186 : sys.solutionOld().get(_dof_indices, &_u_old[0]);
166 464524 : if (_need_u_older)
167 2148 : sys.solutionOlder().get(_dof_indices, &_u_older[0]);
168 :
169 465644 : for (auto tag : _required_vector_tags)
170 2240 : if ((sys.subproblem().vectorTagType(tag) == Moose::VECTOR_TAG_RESIDUAL &&
171 2240 : safe_access_tagged_vectors) ||
172 0 : sys.subproblem().vectorTagType(tag) == Moose::VECTOR_TAG_SOLUTION)
173 1120 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
174 448 : sys.getVector(tag).get(_dof_indices, &_vector_tag_u[tag][0]);
175 :
176 464524 : if (safe_access_tagged_matrices)
177 394710 : for (auto tag : active_coupleable_matrix_tags)
178 448 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
179 224 : for (std::size_t i = 0; i != n; ++i)
180 112 : _matrix_tag_u[tag][i] = sys.getMatrix(tag)(_dof_indices[i], _dof_indices[i]);
181 :
182 464524 : if (_need_u_dot)
183 254012 : (*u_dot).get(_dof_indices, &_u_dot[0]);
184 :
185 464524 : if (_need_u_dotdot)
186 0 : (*u_dotdot).get(_dof_indices, &_u_dotdot[0]);
187 :
188 464524 : if (_need_u_dot_old)
189 0 : (*u_dot_old).get(_dof_indices, &_u_dot_old[0]);
190 :
191 464524 : if (_need_u_dotdot_old)
192 0 : (*u_dotdot_old).get(_dof_indices, &_u_dotdot_old[0]);
193 : }
194 : else
195 : {
196 0 : for (std::size_t i = 0; i != n; ++i)
197 : {
198 0 : const dof_id_type dof_index = _dof_indices[i];
199 0 : std::vector<dof_id_type> one_dof_index(1, dof_index);
200 0 : if (_dof_map.all_semilocal_indices(one_dof_index))
201 : {
202 : libmesh_assert_less(i, _u.size());
203 :
204 0 : current_solution.get(one_dof_index, &_u[i]);
205 0 : if (_need_u_old)
206 0 : sys.solutionOld().get(one_dof_index, &_u_old[i]);
207 0 : if (_need_u_older)
208 0 : sys.solutionOlder().get(one_dof_index, &_u_older[i]);
209 :
210 0 : if (safe_access_tagged_vectors)
211 : {
212 0 : for (auto tag : _required_vector_tags)
213 0 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
214 0 : sys.getVector(tag).get(one_dof_index, &_vector_tag_u[tag][i]);
215 : }
216 :
217 0 : if (safe_access_tagged_matrices)
218 0 : for (auto tag : active_coupleable_matrix_tags)
219 0 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
220 0 : _matrix_tag_u[tag][i] = sys.getMatrix(tag)(dof_index, dof_index);
221 :
222 0 : if (_need_u_dot)
223 0 : (*u_dot).get(one_dof_index, &_u_dot[i]);
224 :
225 0 : if (_need_u_dotdot)
226 0 : (*u_dotdot).get(one_dof_index, &_u_dotdot[i]);
227 :
228 0 : if (_need_u_dot_old)
229 0 : (*u_dot_old).get(one_dof_index, &_u_dot_old[i]);
230 :
231 0 : if (_need_u_dotdot_old)
232 0 : (*u_dotdot_old).get(one_dof_index, &_u_dotdot_old[i]);
233 : }
234 : else
235 : {
236 : #ifdef _GLIBCXX_DEBUG
237 : // Let's make it possible to catch invalid accesses to these
238 : // variables immediately via a thrown exception, if our
239 : // libstdc++ compiler flags allow for that.
240 : _u.resize(i);
241 : if (_need_u_old)
242 : _u_old.resize(i);
243 : if (_need_u_older)
244 : _u_older.resize(i);
245 :
246 : for (auto tag : _required_vector_tags)
247 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
248 : _vector_tag_u[tag].resize(i);
249 :
250 : for (auto tag : active_coupleable_matrix_tags)
251 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
252 : _matrix_tag_u[tag].resize(i);
253 :
254 : if (_need_u_dot)
255 : _u_dot.resize(i);
256 :
257 : if (_need_u_dotdot)
258 : _u_dotdot.resize(i);
259 :
260 : if (_need_u_dot_old)
261 : _u_dot_old.resize(i);
262 :
263 : if (_need_u_dotdot_old)
264 : _u_dotdot_old.resize(i);
265 : #else
266 : // If we can't catch errors at run-time, we can at least
267 : // propagate NaN values rather than invalid values, so that
268 : // users won't trust the result.
269 0 : _u[i] = std::numeric_limits<Real>::quiet_NaN();
270 0 : if (_need_u_old)
271 0 : _u_old[i] = std::numeric_limits<Real>::quiet_NaN();
272 0 : if (_need_u_older)
273 0 : _u_older[i] = std::numeric_limits<Real>::quiet_NaN();
274 :
275 0 : for (auto tag : _required_vector_tags)
276 0 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
277 0 : _vector_tag_u[tag][i] = std::numeric_limits<Real>::quiet_NaN();
278 :
279 0 : for (auto tag : active_coupleable_matrix_tags)
280 0 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
281 0 : _matrix_tag_u[tag][i] = std::numeric_limits<Real>::quiet_NaN();
282 :
283 0 : if (_need_u_dot)
284 0 : _u_dot[i] = std::numeric_limits<Real>::quiet_NaN();
285 :
286 0 : if (_need_u_dotdot)
287 0 : _u_dotdot[i] = std::numeric_limits<Real>::quiet_NaN();
288 :
289 0 : if (_need_u_dot_old)
290 0 : _u_dot_old[i] = std::numeric_limits<Real>::quiet_NaN();
291 :
292 0 : if (_need_u_dotdot_old)
293 0 : _u_dotdot_old[i] = std::numeric_limits<Real>::quiet_NaN();
294 : #endif
295 : }
296 0 : }
297 : }
298 :
299 : // Automatic differentiation
300 464524 : if (_need_ad)
301 18723 : computeAD(/*nodal_ordering=*/false);
302 : }
303 :
304 : void
305 19316 : MooseVariableScalar::computeAD(bool)
306 : {
307 19316 : auto n_dofs = _dof_indices.size();
308 : const bool do_derivatives =
309 19316 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
310 :
311 19316 : if (_need_ad_u)
312 : {
313 19316 : _ad_u.resize(n_dofs);
314 39678 : for (MooseIndex(n_dofs) i = 0; i < n_dofs; ++i)
315 : {
316 20362 : _ad_u[i] = _u[i];
317 20362 : if (do_derivatives)
318 9074 : Moose::derivInsert(_ad_u[i].derivatives(), _dof_indices[i], 1.);
319 : }
320 : }
321 :
322 19316 : if (_need_ad_u_dot)
323 : {
324 1839 : _ad_u_dot.resize(n_dofs);
325 3678 : for (MooseIndex(n_dofs) i = 0; i < n_dofs; ++i)
326 : {
327 1839 : _ad_u_dot[i] = _u_dot[i];
328 1839 : if (do_derivatives)
329 1047 : Moose::derivInsert(_ad_u_dot[i].derivatives(), _dof_indices[i], _du_dot_du[i]);
330 : }
331 : }
332 19316 : }
333 :
334 : void
335 19302 : MooseVariableScalar::setValue(unsigned int i, Number value)
336 : {
337 : // In debug modes, we might have set a "trap" to catch reads of
338 : // uninitialized values, but this trap shouldn't prevent setting
339 : // values.
340 : #ifdef DEBUG
341 : if (i >= _u.size())
342 : {
343 : libmesh_assert_less(i, _dof_indices.size());
344 : _u.resize(i + 1);
345 : }
346 : #endif
347 19302 : _u[i] = value; // update variable value
348 19302 : }
349 :
350 : void
351 206 : MooseVariableScalar::setValues(Number value)
352 : {
353 206 : unsigned int n = _dof_indices.size();
354 : // In debug modes, we might have set a "trap" to catch reads of
355 : // uninitialized values, but this trap shouldn't prevent setting
356 : // values.
357 : #ifdef DEBUG
358 : _u.resize(n);
359 : #endif
360 412 : for (unsigned int i = 0; i < n; i++)
361 206 : _u[i] = value;
362 206 : }
363 :
364 : void
365 29384 : MooseVariableScalar::insert(NumericVector<Number> & soln)
366 : {
367 : // We may have redundantly computed this value on many different
368 : // processors, but only the processor which actually owns it should
369 : // be saving it to the solution vector, to avoid O(N_scalar_vars)
370 : // unnecessary communication.
371 :
372 29384 : const dof_id_type first_dof = _dof_map.first_dof();
373 29384 : const dof_id_type end_dof = _dof_map.end_dof();
374 29384 : if (_dof_indices.size() > 0 && first_dof <= _dof_indices[0] && _dof_indices[0] < end_dof)
375 25459 : soln.insert(&_u[0], _dof_indices);
376 29384 : }
377 :
378 : const ADVariableValue &
379 577 : MooseVariableScalar::adSln() const
380 : {
381 577 : _need_ad = _need_ad_u = true;
382 577 : return _ad_u;
383 : }
384 :
385 : const VariableValue &
386 39 : MooseVariableScalar::slnOld() const
387 : {
388 39 : _need_u_old = true;
389 39 : return _u_old;
390 : }
391 :
392 : const VariableValue &
393 26 : MooseVariableScalar::slnOlder() const
394 : {
395 26 : _need_u_older = true;
396 26 : return _u_older;
397 : }
398 :
399 : const VariableValue &
400 28 : MooseVariableScalar::vectorTagSln(TagID tag) const
401 : {
402 28 : _need_vector_tag_u[tag] = true;
403 28 : return _vector_tag_u[tag];
404 : }
405 :
406 : const VariableValue &
407 14 : MooseVariableScalar::matrixTagSln(TagID tag) const
408 : {
409 14 : _need_matrix_tag_u[tag] = true;
410 14 : return _matrix_tag_u[tag];
411 : }
412 :
413 : const ADVariableValue &
414 58 : MooseVariableScalar::adUDot() const
415 : {
416 58 : if (_sys.solutionUDot())
417 : {
418 58 : _need_ad = _need_ad_u_dot = _need_u_dot = true;
419 58 : return _ad_u_dot;
420 : }
421 : else
422 0 : mooseError("MooseVariableScalar: Time derivative of solution (`u_dot`) is not stored. Please "
423 : "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
424 : }
425 :
426 : const VariableValue &
427 466 : MooseVariableScalar::uDot() const
428 : {
429 466 : if (_sys.solutionUDot())
430 : {
431 466 : _need_u_dot = true;
432 466 : return _u_dot;
433 : }
434 : else
435 0 : mooseError("MooseVariableScalar: Time derivative of solution (`u_dot`) is not stored. Please "
436 : "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
437 : }
438 :
439 : const VariableValue &
440 0 : MooseVariableScalar::uDotDot() const
441 : {
442 0 : if (_sys.solutionUDotDot())
443 : {
444 0 : _need_u_dotdot = true;
445 0 : return _u_dotdot;
446 : }
447 : else
448 0 : mooseError("MooseVariableScalar: Second time derivative of solution (`u_dotdot`) is not "
449 : "stored. Please set uDotDotRequested() to true in FEProblemBase before requesting "
450 : "`u_dotdot`.");
451 : }
452 :
453 : const VariableValue &
454 0 : MooseVariableScalar::uDotOld() const
455 : {
456 0 : if (_sys.solutionUDotOld())
457 : {
458 0 : _need_u_dot_old = true;
459 0 : return _u_dot_old;
460 : }
461 : else
462 0 : mooseError("MooseVariableScalar: Old time derivative of solution (`u_dot_old`) is not "
463 : "stored. Please set uDotOldRequested() to true in FEProblemBase before requesting "
464 : "`u_dot_old`.");
465 : }
466 :
467 : const VariableValue &
468 0 : MooseVariableScalar::uDotDotOld() const
469 : {
470 0 : if (_sys.solutionUDotDotOld())
471 : {
472 0 : _need_u_dotdot_old = true;
473 0 : return _u_dotdot_old;
474 : }
475 : else
476 0 : mooseError("MooseVariableScalar: Old second time derivative of solution (`u_dotdot_old`) is "
477 : "not stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
478 : "requesting `u_dotdot_old`.");
479 : }
480 :
481 : const VariableValue &
482 466 : MooseVariableScalar::duDotDu() const
483 : {
484 466 : _need_du_dot_du = true;
485 466 : return _du_dot_du;
486 : }
487 :
488 : const VariableValue &
489 0 : MooseVariableScalar::duDotDotDu() const
490 : {
491 0 : _need_du_dotdot_du = true;
492 0 : return _du_dotdot_du;
493 : }
494 :
495 : unsigned int
496 2804 : MooseVariableScalar::oldestSolutionStateRequested() const
497 : {
498 2804 : if (_need_u_older)
499 12 : return 2;
500 2792 : if (_need_u_old)
501 24 : return 1;
502 2768 : return 0;
503 : }
504 :
505 : MooseVariableScalar::ValueType
506 0 : MooseVariableScalar::evaluate(const ElemArg & /*elem_arg*/, const Moose::StateArg & state) const
507 : {
508 0 : return evaluate(state);
509 : }
510 :
511 : MooseVariableScalar::ValueType
512 0 : MooseVariableScalar::evaluate(const FaceArg & /*face*/, const Moose::StateArg & state) const
513 : {
514 0 : return evaluate(state);
515 : }
516 :
517 : MooseVariableScalar::ValueType
518 150 : MooseVariableScalar::evaluate(const ElemQpArg & /*elem_qp*/, const Moose::StateArg & state) const
519 : {
520 150 : return evaluate(state);
521 : }
522 :
523 : MooseVariableScalar::ValueType
524 0 : MooseVariableScalar::evaluate(const ElemSideQpArg & /*elem_side_qp*/,
525 : const Moose::StateArg & state) const
526 : {
527 0 : return evaluate(state);
528 : }
529 :
530 : MooseVariableScalar::ValueType
531 0 : MooseVariableScalar::evaluate(const ElemPointArg & /*elem_point_arg*/,
532 : const Moose::StateArg & state) const
533 : {
534 0 : return evaluate(state);
535 : }
536 :
537 : MooseVariableScalar::ValueType
538 0 : MooseVariableScalar::evaluate(const NodeArg & /*node_arg*/, const Moose::StateArg & state) const
539 : {
540 0 : return evaluate(state);
541 : }
542 :
543 : MooseVariableScalar::ValueType
544 150 : MooseVariableScalar::evaluate(const Moose::StateArg & state) const
545 : {
546 : mooseAssert(_dof_indices.size() == 1,
547 : "MooseVariableScalar::evaluate may only be used for first-order scalar variables");
548 :
549 150 : const auto & solution = getSolution(state);
550 150 : ADReal value = solution(_dof_indices[0]);
551 150 : if (state.state == 0 && Moose::doDerivatives(_subproblem, _sys))
552 50 : Moose::derivInsert(value.derivatives(), _dof_indices[0], 1.0);
553 :
554 150 : return value;
555 0 : }
556 :
557 : MooseVariableScalar::GradientType
558 0 : MooseVariableScalar::evaluateGradient(const ElemArg & /*elem_arg*/,
559 : const Moose::StateArg & /*state*/) const
560 : {
561 0 : return 0;
562 : }
563 :
564 : MooseVariableScalar::GradientType
565 0 : MooseVariableScalar::evaluateGradient(const FaceArg & /*face*/,
566 : const Moose::StateArg & /*state*/) const
567 : {
568 0 : return 0;
569 : }
570 :
571 : MooseVariableScalar::GradientType
572 0 : MooseVariableScalar::evaluateGradient(const ElemQpArg & /*elem_qp*/,
573 : const Moose::StateArg & /*state*/) const
574 : {
575 0 : return 0;
576 : }
577 :
578 : MooseVariableScalar::GradientType
579 0 : MooseVariableScalar::evaluateGradient(const ElemSideQpArg & /*elem_side_qp*/,
580 : const Moose::StateArg & /*state*/) const
581 : {
582 0 : return 0;
583 : }
584 :
585 : MooseVariableScalar::GradientType
586 0 : MooseVariableScalar::evaluateGradient(const ElemPointArg & /*elem_point_arg*/,
587 : const Moose::StateArg & /*state*/) const
588 : {
589 0 : return 0;
590 : }
591 :
592 : MooseVariableScalar::GradientType
593 0 : MooseVariableScalar::evaluateGradient(const NodeArg & /*node_arg*/,
594 : const Moose::StateArg & /*state*/) const
595 : {
596 0 : return 0;
597 : }
|
