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