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 17711 : MooseVariableScalar::validParams()
26 : {
27 17711 : auto params = MooseVariableBase::validParams();
28 17711 : params.addClassDescription("Moose wrapper class around scalar variables");
29 17711 : params.set<MooseEnum>("family") = "SCALAR";
30 17711 : return params;
31 0 : }
32 :
33 3290 : MooseVariableScalar::MooseVariableScalar(const InputParameters & parameters)
34 : : MooseVariableBase(parameters),
35 3290 : _need_u_dot(false),
36 3290 : _need_u_dotdot(false),
37 3290 : _need_u_dot_old(false),
38 3290 : _need_u_dotdot_old(false),
39 3290 : _need_du_dot_du(false),
40 3290 : _need_du_dotdot_du(false),
41 3290 : _need_u_old(false),
42 3290 : _need_u_older(false),
43 3290 : _need_ad(false),
44 3290 : _need_ad_u(false),
45 3290 : _need_ad_u_dot(false)
46 : {
47 3290 : auto num_vector_tags = _sys.subproblem().numVectorTags();
48 :
49 3290 : _vector_tag_u.resize(num_vector_tags);
50 3290 : _need_vector_tag_u.resize(num_vector_tags);
51 :
52 3290 : auto num_matrix_tags = _sys.subproblem().numMatrixTags();
53 :
54 3290 : _matrix_tag_u.resize(num_matrix_tags);
55 3290 : _need_matrix_tag_u.resize(num_matrix_tags);
56 3290 : }
57 :
58 6218 : MooseVariableScalar::~MooseVariableScalar()
59 : {
60 3109 : _u.release();
61 3109 : _u_old.release();
62 3109 : _u_older.release();
63 :
64 3109 : _u_dot.release();
65 3109 : _u_dotdot.release();
66 3109 : _u_dot_old.release();
67 3109 : _u_dotdot_old.release();
68 3109 : _du_dot_du.release();
69 3109 : _du_dotdot_du.release();
70 :
71 21312 : for (auto & _tag_u : _vector_tag_u)
72 18203 : _tag_u.release();
73 :
74 3109 : _vector_tag_u.clear();
75 :
76 9549 : for (auto & _tag_u : _matrix_tag_u)
77 6440 : _tag_u.release();
78 :
79 3109 : _matrix_tag_u.clear();
80 6218 : }
81 :
82 : void
83 539212 : MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/)
84 : {
85 539212 : if (reinit_for_derivative_reordering)
86 : {
87 : // We've already calculated everything in an earlier reinit. All we have to do is re-sort the
88 : // derivative vector for AD calculations
89 18703 : if (_need_ad)
90 602 : computeAD(/*nodal_ordering=*/true);
91 18703 : return;
92 : }
93 :
94 : // We want to make sure that we're not allocating data with any of the
95 : // accessors that follow, but _sys is non-const
96 520509 : const SystemBase & sys = _sys;
97 :
98 520509 : const NumericVector<Real> & current_solution = *sys.currentSolution();
99 520509 : const NumericVector<Real> * u_dot = sys.solutionUDot();
100 520509 : const NumericVector<Real> * u_dotdot = sys.solutionUDotDot();
101 520509 : const NumericVector<Real> * u_dot_old = sys.solutionUDotOld();
102 520509 : const NumericVector<Real> * u_dotdot_old = sys.solutionUDotDotOld();
103 520509 : const Real & du_dot_du = sys.duDotDu(number());
104 520509 : const Real & du_dotdot_du = sys.duDotDotDu();
105 520509 : auto safe_access_tagged_vectors = sys.subproblem().safeAccessTaggedVectors();
106 520509 : auto safe_access_tagged_matrices = sys.subproblem().safeAccessTaggedMatrices();
107 : auto & active_coupleable_matrix_tags =
108 520509 : sys.subproblem().getActiveScalarVariableCoupleableMatrixTags(_tid);
109 :
110 520509 : _dof_map.SCALAR_dof_indices(_dof_indices, _var_num);
111 :
112 520509 : auto n = _dof_indices.size();
113 520509 : _u.resize(n);
114 520509 : if (_need_u_old)
115 8094 : _u_old.resize(n);
116 520509 : if (_need_u_older)
117 2108 : _u_older.resize(n);
118 :
119 4527258 : for (auto & _tag_u : _vector_tag_u)
120 4006749 : _tag_u.resize(n);
121 :
122 1576375 : for (auto & _tag_u : _matrix_tag_u)
123 1055866 : _tag_u.resize(n);
124 :
125 520509 : _du_dot_du.clear();
126 520509 : _du_dot_du.resize(n, du_dot_du);
127 :
128 520509 : if (_need_u_dot)
129 286133 : _u_dot.resize(n);
130 :
131 520509 : if (_need_u_dotdot)
132 0 : _u_dotdot.resize(n);
133 :
134 520509 : if (_need_u_dot_old)
135 0 : _u_dot_old.resize(n);
136 :
137 520509 : if (_need_u_dotdot_old)
138 0 : _u_dotdot_old.resize(n);
139 :
140 520509 : if (_need_du_dot_du)
141 : {
142 284403 : _du_dot_du.clear();
143 284403 : _du_dot_du.resize(n, du_dot_du);
144 : }
145 :
146 520509 : if (_need_du_dotdot_du)
147 : {
148 0 : _du_dotdot_du.clear();
149 0 : _du_dotdot_du.resize(n, du_dotdot_du);
150 : }
151 :
152 : // If we have an empty partition, or if we have a partition which
153 : // does not include any of the subdomains of a subdomain-restricted
154 : // variable, then we do not have access to that variable! Hopefully
155 : // we won't need the indices we lack.
156 520509 : if (_dof_map.all_semilocal_indices(_dof_indices))
157 : {
158 520509 : current_solution.get(_dof_indices, &_u[0]);
159 520509 : if (_need_u_old)
160 8094 : sys.solutionOld().get(_dof_indices, &_u_old[0]);
161 520509 : if (_need_u_older)
162 2108 : sys.solutionOlder().get(_dof_indices, &_u_older[0]);
163 :
164 522009 : for (auto tag : _required_vector_tags)
165 3000 : if ((sys.subproblem().vectorTagType(tag) == Moose::VECTOR_TAG_RESIDUAL &&
166 3000 : safe_access_tagged_vectors) ||
167 0 : sys.subproblem().vectorTagType(tag) == Moose::VECTOR_TAG_SOLUTION)
168 1500 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
169 600 : sys.getVector(tag).get(_dof_indices, &_vector_tag_u[tag][0]);
170 :
171 520509 : if (safe_access_tagged_matrices)
172 : {
173 441559 : for (auto tag : active_coupleable_matrix_tags)
174 600 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
175 300 : for (std::size_t i = 0; i != n; ++i)
176 : {
177 150 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
178 150 : _matrix_tag_u[tag][i] = sys.getMatrix(tag)(_dof_indices[i], _dof_indices[i]);
179 150 : }
180 : }
181 :
182 520509 : if (_need_u_dot)
183 286133 : (*u_dot).get(_dof_indices, &_u_dot[0]);
184 :
185 520509 : if (_need_u_dotdot)
186 0 : (*u_dotdot).get(_dof_indices, &_u_dotdot[0]);
187 :
188 520509 : if (_need_u_dot_old)
189 0 : (*u_dot_old).get(_dof_indices, &_u_dot_old[0]);
190 :
191 520509 : 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 : {
219 0 : for (auto tag : active_coupleable_matrix_tags)
220 0 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
221 : {
222 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
223 0 : _matrix_tag_u[tag][i] = sys.getMatrix(tag)(dof_index, dof_index);
224 0 : }
225 : }
226 :
227 0 : if (_need_u_dot)
228 0 : (*u_dot).get(one_dof_index, &_u_dot[i]);
229 :
230 0 : if (_need_u_dotdot)
231 0 : (*u_dotdot).get(one_dof_index, &_u_dotdot[i]);
232 :
233 0 : if (_need_u_dot_old)
234 0 : (*u_dot_old).get(one_dof_index, &_u_dot_old[i]);
235 :
236 0 : if (_need_u_dotdot_old)
237 0 : (*u_dotdot_old).get(one_dof_index, &_u_dotdot_old[i]);
238 : }
239 : else
240 : {
241 : #ifdef _GLIBCXX_DEBUG
242 : // Let's make it possible to catch invalid accesses to these
243 : // variables immediately via a thrown exception, if our
244 : // libstdc++ compiler flags allow for that.
245 : _u.resize(i);
246 : if (_need_u_old)
247 : _u_old.resize(i);
248 : if (_need_u_older)
249 : _u_older.resize(i);
250 :
251 : for (auto tag : _required_vector_tags)
252 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
253 : _vector_tag_u[tag].resize(i);
254 :
255 : for (auto tag : active_coupleable_matrix_tags)
256 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
257 : _matrix_tag_u[tag].resize(i);
258 :
259 : if (_need_u_dot)
260 : _u_dot.resize(i);
261 :
262 : if (_need_u_dotdot)
263 : _u_dotdot.resize(i);
264 :
265 : if (_need_u_dot_old)
266 : _u_dot_old.resize(i);
267 :
268 : if (_need_u_dotdot_old)
269 : _u_dotdot_old.resize(i);
270 : #else
271 : // If we can't catch errors at run-time, we can at least
272 : // propagate NaN values rather than invalid values, so that
273 : // users won't trust the result.
274 0 : _u[i] = std::numeric_limits<Real>::quiet_NaN();
275 0 : if (_need_u_old)
276 0 : _u_old[i] = std::numeric_limits<Real>::quiet_NaN();
277 0 : if (_need_u_older)
278 0 : _u_older[i] = std::numeric_limits<Real>::quiet_NaN();
279 :
280 0 : for (auto tag : _required_vector_tags)
281 0 : if (sys.hasVector(tag) && _need_vector_tag_u[tag])
282 0 : _vector_tag_u[tag][i] = std::numeric_limits<Real>::quiet_NaN();
283 :
284 0 : for (auto tag : active_coupleable_matrix_tags)
285 0 : if (sys.hasMatrix(tag) && sys.getMatrix(tag).closed() && _need_matrix_tag_u[tag])
286 0 : _matrix_tag_u[tag][i] = std::numeric_limits<Real>::quiet_NaN();
287 :
288 0 : if (_need_u_dot)
289 0 : _u_dot[i] = std::numeric_limits<Real>::quiet_NaN();
290 :
291 0 : if (_need_u_dotdot)
292 0 : _u_dotdot[i] = std::numeric_limits<Real>::quiet_NaN();
293 :
294 0 : if (_need_u_dot_old)
295 0 : _u_dot_old[i] = std::numeric_limits<Real>::quiet_NaN();
296 :
297 0 : if (_need_u_dotdot_old)
298 0 : _u_dotdot_old[i] = std::numeric_limits<Real>::quiet_NaN();
299 : #endif
300 : }
301 0 : }
302 : }
303 :
304 : // Automatic differentiation
305 520509 : if (_need_ad)
306 19093 : computeAD(/*nodal_ordering=*/false);
307 : }
308 :
309 : void
310 19695 : MooseVariableScalar::computeAD(bool)
311 : {
312 19695 : auto n_dofs = _dof_indices.size();
313 : const bool do_derivatives =
314 19695 : ADReal::do_derivatives && _sys.number() == _subproblem.currentNlSysNum();
315 :
316 19695 : if (_need_ad_u)
317 : {
318 19695 : _ad_u.resize(n_dofs);
319 40426 : for (MooseIndex(n_dofs) i = 0; i < n_dofs; ++i)
320 : {
321 20731 : _ad_u[i] = _u[i];
322 20731 : if (do_derivatives)
323 9189 : Moose::derivInsert(_ad_u[i].derivatives(), _dof_indices[i], 1.);
324 : }
325 : }
326 :
327 19695 : if (_need_ad_u_dot)
328 : {
329 1792 : _ad_u_dot.resize(n_dofs);
330 3584 : for (MooseIndex(n_dofs) i = 0; i < n_dofs; ++i)
331 : {
332 1792 : _ad_u_dot[i] = _u_dot[i];
333 1792 : if (do_derivatives)
334 924 : Moose::derivInsert(_ad_u_dot[i].derivatives(), _dof_indices[i], _du_dot_du[i]);
335 : }
336 : }
337 19695 : }
338 :
339 : void
340 20542 : MooseVariableScalar::setValue(unsigned int i, Number value)
341 : {
342 : // In debug modes, we might have set a "trap" to catch reads of
343 : // uninitialized values, but this trap shouldn't prevent setting
344 : // values.
345 : #ifdef DEBUG
346 : if (i >= _u.size())
347 : {
348 : libmesh_assert_less(i, _dof_indices.size());
349 : _u.resize(i + 1);
350 : }
351 : #endif
352 20542 : _u[i] = value; // update variable value
353 20542 : }
354 :
355 : void
356 1416 : MooseVariableScalar::setValues(Number value)
357 : {
358 1416 : unsigned int n = _dof_indices.size();
359 : // In debug modes, we might have set a "trap" to catch reads of
360 : // uninitialized values, but this trap shouldn't prevent setting
361 : // values.
362 : #ifdef DEBUG
363 : _u.resize(n);
364 : #endif
365 2832 : for (unsigned int i = 0; i < n; i++)
366 1416 : _u[i] = value;
367 1416 : }
368 :
369 : void
370 35625 : MooseVariableScalar::insert(NumericVector<Number> & soln)
371 : {
372 : // We may have redundantly computed this value on many different
373 : // processors, but only the processor which actually owns it should
374 : // be saving it to the solution vector, to avoid O(N_scalar_vars)
375 : // unnecessary communication.
376 :
377 35625 : const dof_id_type first_dof = _dof_map.first_dof();
378 35625 : const dof_id_type end_dof = _dof_map.end_dof();
379 35625 : if (_dof_indices.size() > 0 && first_dof <= _dof_indices[0] && _dof_indices[0] < end_dof)
380 31737 : soln.insert(&_u[0], _dof_indices);
381 35625 : }
382 :
383 : const ADVariableValue &
384 592 : MooseVariableScalar::adSln() const
385 : {
386 592 : _need_ad = _need_ad_u = true;
387 592 : return _ad_u;
388 : }
389 :
390 : const VariableValue &
391 39 : MooseVariableScalar::slnOld() const
392 : {
393 39 : _need_u_old = true;
394 39 : return _u_old;
395 : }
396 :
397 : const VariableValue &
398 26 : MooseVariableScalar::slnOlder() const
399 : {
400 26 : _need_u_older = true;
401 26 : return _u_older;
402 : }
403 :
404 : const VariableValue &
405 36 : MooseVariableScalar::vectorTagSln(TagID tag) const
406 : {
407 36 : _need_vector_tag_u[tag] = true;
408 36 : return _vector_tag_u[tag];
409 : }
410 :
411 : const VariableValue &
412 18 : MooseVariableScalar::matrixTagSln(TagID tag) const
413 : {
414 18 : _need_matrix_tag_u[tag] = true;
415 18 : return _matrix_tag_u[tag];
416 : }
417 :
418 : const ADVariableValue &
419 56 : MooseVariableScalar::adUDot() const
420 : {
421 56 : if (_sys.solutionUDot())
422 : {
423 56 : _need_ad = _need_ad_u_dot = _need_u_dot = true;
424 56 : return _ad_u_dot;
425 : }
426 : else
427 0 : mooseError("MooseVariableScalar: Time derivative of solution (`u_dot`) is not stored. Please "
428 : "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
429 : }
430 :
431 : const VariableValue &
432 515 : MooseVariableScalar::uDot() const
433 : {
434 515 : if (_sys.solutionUDot())
435 : {
436 515 : _need_u_dot = true;
437 515 : return _u_dot;
438 : }
439 : else
440 0 : mooseError("MooseVariableScalar: Time derivative of solution (`u_dot`) is not stored. Please "
441 : "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
442 : }
443 :
444 : const VariableValue &
445 0 : MooseVariableScalar::uDotDot() const
446 : {
447 0 : if (_sys.solutionUDotDot())
448 : {
449 0 : _need_u_dotdot = true;
450 0 : return _u_dotdot;
451 : }
452 : else
453 0 : mooseError("MooseVariableScalar: Second time derivative of solution (`u_dotdot`) is not "
454 : "stored. Please set uDotDotRequested() to true in FEProblemBase before requesting "
455 : "`u_dotdot`.");
456 : }
457 :
458 : const VariableValue &
459 0 : MooseVariableScalar::uDotOld() const
460 : {
461 0 : if (_sys.solutionUDotOld())
462 : {
463 0 : _need_u_dot_old = true;
464 0 : return _u_dot_old;
465 : }
466 : else
467 0 : mooseError("MooseVariableScalar: Old time derivative of solution (`u_dot_old`) is not "
468 : "stored. Please set uDotOldRequested() to true in FEProblemBase before requesting "
469 : "`u_dot_old`.");
470 : }
471 :
472 : const VariableValue &
473 0 : MooseVariableScalar::uDotDotOld() const
474 : {
475 0 : if (_sys.solutionUDotDotOld())
476 : {
477 0 : _need_u_dotdot_old = true;
478 0 : return _u_dotdot_old;
479 : }
480 : else
481 0 : mooseError("MooseVariableScalar: Old second time derivative of solution (`u_dotdot_old`) is "
482 : "not stored. Please set uDotDotOldRequested() to true in FEProblemBase before "
483 : "requesting `u_dotdot_old`.");
484 : }
485 :
486 : const VariableValue &
487 515 : MooseVariableScalar::duDotDu() const
488 : {
489 515 : _need_du_dot_du = true;
490 515 : return _du_dot_du;
491 : }
492 :
493 : const VariableValue &
494 0 : MooseVariableScalar::duDotDotDu() const
495 : {
496 0 : _need_du_dotdot_du = true;
497 0 : return _du_dotdot_du;
498 : }
499 :
500 : unsigned int
501 2954 : MooseVariableScalar::oldestSolutionStateRequested() const
502 : {
503 2954 : if (_need_u_older)
504 12 : return 2;
505 2942 : if (_need_u_old)
506 24 : return 1;
507 2918 : return 0;
508 : }
|
