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 "MooseVariableFV.h"
11 : #include "TimeIntegrator.h"
12 : #include "NonlinearSystemBase.h"
13 : #include "DisplacedSystem.h"
14 : #include "SystemBase.h"
15 : #include "SubProblem.h"
16 : #include "Assembly.h"
17 : #include "MathFVUtils.h"
18 : #include "FVUtils.h"
19 : #include "FVFluxBC.h"
20 : #include "FVDirichletBCBase.h"
21 : #include "GreenGaussGradient.h"
22 :
23 : #include "libmesh/numeric_vector.h"
24 :
25 : #include <climits>
26 : #include <typeinfo>
27 :
28 : using namespace Moose;
29 :
30 : registerMooseObject("MooseApp", MooseVariableFVReal);
31 :
32 : template <typename OutputType>
33 : InputParameters
34 12022 : MooseVariableFV<OutputType>::validParams()
35 : {
36 12022 : InputParameters params = MooseVariableField<OutputType>::validParams();
37 24044 : params.set<bool>("fv") = true;
38 48088 : params.set<MooseEnum>("family") = "MONOMIAL";
39 48088 : params.set<MooseEnum>("order") = "CONSTANT";
40 36066 : params.template addParam<bool>(
41 : "two_term_boundary_expansion",
42 24044 : true,
43 : "Whether to use a two-term Taylor expansion to calculate boundary face values. "
44 : "If the two-term expansion is used, then the boundary face value depends on the "
45 : "adjoining cell center gradient, which itself depends on the boundary face value. "
46 : "Consequently an implicit solve is used to simultaneously solve for the adjoining cell "
47 : "center gradient and boundary face value(s).");
48 48088 : MooseEnum face_interp_method("average skewness-corrected", "average");
49 48088 : params.template addParam<MooseEnum>("face_interp_method",
50 : face_interp_method,
51 : "Switch that can select between face interpolation methods.");
52 24044 : params.template addParam<bool>(
53 24044 : "cache_cell_gradients", true, "Whether to cache cell gradients or re-compute them.");
54 :
55 : // Depending on the face interpolation we might have to do more than one layer ghosting.
56 36066 : params.addRelationshipManager(
57 : "ElementSideNeighborLayers",
58 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
59 : Moose::RelationshipManagerType::COUPLING,
60 0 : [](const InputParameters & obj_params, InputParameters & rm_params)
61 : {
62 10145 : unsigned short layers = 1;
63 10145 : if (obj_params.get<MooseEnum>("face_interp_method") == "skewness-corrected")
64 180 : layers = 2;
65 :
66 30435 : rm_params.set<unsigned short>("layers") = layers;
67 : });
68 24044 : return params;
69 12022 : }
70 :
71 : template <typename OutputType>
72 5421 : MooseVariableFV<OutputType>::MooseVariableFV(const InputParameters & parameters)
73 : : MooseVariableField<OutputType>(parameters),
74 5421 : _solution(this->_sys.currentSolution()),
75 5421 : _phi(this->_assembly.template fePhi<OutputShape>(FEType(CONSTANT, MONOMIAL))),
76 5421 : _grad_phi(this->_assembly.template feGradPhi<OutputShape>(FEType(CONSTANT, MONOMIAL))),
77 5421 : _phi_face(this->_assembly.template fePhiFace<OutputShape>(FEType(CONSTANT, MONOMIAL))),
78 5421 : _grad_phi_face(this->_assembly.template feGradPhiFace<OutputShape>(FEType(CONSTANT, MONOMIAL))),
79 5421 : _phi_face_neighbor(
80 5421 : this->_assembly.template fePhiFaceNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
81 5421 : _grad_phi_face_neighbor(
82 5421 : this->_assembly.template feGradPhiFaceNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
83 5421 : _phi_neighbor(this->_assembly.template fePhiNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
84 5421 : _grad_phi_neighbor(
85 5421 : this->_assembly.template feGradPhiNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
86 5421 : _prev_elem(nullptr),
87 10842 : _two_term_boundary_expansion(this->isParamValid("two_term_boundary_expansion")
88 16878 : ? this->template getParam<bool>("two_term_boundary_expansion")
89 : : true),
90 10842 : _cache_cell_gradients(this->isParamValid("cache_cell_gradients")
91 16878 : ? this->template getParam<bool>("cache_cell_gradients")
92 16263 : : true)
93 : {
94 5421 : _element_data = std::make_unique<MooseVariableDataFV<OutputType>>(
95 5421 : *this, _sys, _tid, Moose::ElementType::Element, this->_assembly.elem());
96 5421 : _neighbor_data = std::make_unique<MooseVariableDataFV<OutputType>>(
97 5421 : *this, _sys, _tid, Moose::ElementType::Neighbor, this->_assembly.neighbor());
98 :
99 16263 : if (this->isParamValid("face_interp_method"))
100 : {
101 7638 : const auto & interp_method = this->template getParam<MooseEnum>("face_interp_method");
102 3819 : if (interp_method == "average")
103 3759 : _face_interp_method = Moose::FV::InterpMethod::Average;
104 60 : else if (interp_method == "skewness-corrected")
105 60 : _face_interp_method = Moose::FV::InterpMethod::SkewCorrectedAverage;
106 : }
107 : else
108 1602 : _face_interp_method = Moose::FV::InterpMethod::Average;
109 13059 : }
110 :
111 : template <typename OutputType>
112 : void
113 10572276 : MooseVariableFV<OutputType>::clearDofIndices()
114 : {
115 10572276 : _element_data->clearDofIndices();
116 10572276 : }
117 :
118 : template <typename OutputType>
119 : typename MooseVariableFV<OutputType>::DofValue
120 0 : MooseVariableFV<OutputType>::getElementalValue(const Elem * elem, unsigned int idx) const
121 : {
122 0 : return _element_data->getElementalValue(elem, Moose::Current, idx);
123 : }
124 :
125 : template <typename OutputType>
126 : typename MooseVariableFV<OutputType>::DofValue
127 0 : MooseVariableFV<OutputType>::getElementalValueOld(const Elem * elem, unsigned int idx) const
128 : {
129 0 : return _element_data->getElementalValue(elem, Moose::Old, idx);
130 : }
131 :
132 : template <typename OutputType>
133 : typename MooseVariableFV<OutputType>::DofValue
134 0 : MooseVariableFV<OutputType>::getElementalValueOlder(const Elem * elem, unsigned int idx) const
135 : {
136 0 : return _element_data->getElementalValue(elem, Moose::Older, idx);
137 : }
138 :
139 : template <typename OutputType>
140 : void
141 318608 : MooseVariableFV<OutputType>::insert(NumericVector<Number> & residual)
142 : {
143 318608 : _element_data->insert(residual);
144 318608 : }
145 :
146 : template <typename OutputType>
147 : void
148 0 : MooseVariableFV<OutputType>::insertLower(NumericVector<Number> &)
149 : {
150 0 : lowerDError();
151 : }
152 :
153 : template <typename OutputType>
154 : void
155 36136 : MooseVariableFV<OutputType>::add(NumericVector<Number> & residual)
156 : {
157 36136 : _element_data->add(residual);
158 36136 : }
159 :
160 : template <typename OutputType>
161 : const typename MooseVariableFV<OutputType>::DofValues &
162 65 : MooseVariableFV<OutputType>::dofValues() const
163 : {
164 65 : return _element_data->dofValues();
165 : }
166 :
167 : template <typename OutputType>
168 : const typename MooseVariableFV<OutputType>::DofValues &
169 0 : MooseVariableFV<OutputType>::dofValuesOld() const
170 : {
171 0 : return _element_data->dofValuesOld();
172 : }
173 :
174 : template <typename OutputType>
175 : const typename MooseVariableFV<OutputType>::DofValues &
176 0 : MooseVariableFV<OutputType>::dofValuesOlder() const
177 : {
178 0 : return _element_data->dofValuesOlder();
179 : }
180 :
181 : template <typename OutputType>
182 : const typename MooseVariableFV<OutputType>::DofValues &
183 0 : MooseVariableFV<OutputType>::dofValuesPreviousNL() const
184 : {
185 0 : return _element_data->dofValuesPreviousNL();
186 : }
187 :
188 : template <typename OutputType>
189 : const typename MooseVariableFV<OutputType>::DofValues &
190 0 : MooseVariableFV<OutputType>::dofValuesNeighbor() const
191 : {
192 0 : return _neighbor_data->dofValues();
193 : }
194 :
195 : template <typename OutputType>
196 : const typename MooseVariableFV<OutputType>::DofValues &
197 0 : MooseVariableFV<OutputType>::dofValuesOldNeighbor() const
198 : {
199 0 : return _neighbor_data->dofValuesOld();
200 : }
201 :
202 : template <typename OutputType>
203 : const typename MooseVariableFV<OutputType>::DofValues &
204 0 : MooseVariableFV<OutputType>::dofValuesOlderNeighbor() const
205 : {
206 0 : return _neighbor_data->dofValuesOlder();
207 : }
208 :
209 : template <typename OutputType>
210 : const typename MooseVariableFV<OutputType>::DofValues &
211 0 : MooseVariableFV<OutputType>::dofValuesPreviousNLNeighbor() const
212 : {
213 0 : return _neighbor_data->dofValuesPreviousNL();
214 : }
215 :
216 : template <typename OutputType>
217 : const typename MooseVariableFV<OutputType>::DofValues &
218 0 : MooseVariableFV<OutputType>::dofValuesDot() const
219 : {
220 0 : return _element_data->dofValuesDot();
221 : }
222 :
223 : template <typename OutputType>
224 : const typename MooseVariableFV<OutputType>::DofValues &
225 0 : MooseVariableFV<OutputType>::dofValuesDotDot() const
226 : {
227 0 : return _element_data->dofValuesDotDot();
228 : }
229 :
230 : template <typename OutputType>
231 : const typename MooseVariableFV<OutputType>::DofValues &
232 0 : MooseVariableFV<OutputType>::dofValuesDotOld() const
233 : {
234 0 : return _element_data->dofValuesDotOld();
235 : }
236 :
237 : template <typename OutputType>
238 : const typename MooseVariableFV<OutputType>::DofValues &
239 0 : MooseVariableFV<OutputType>::dofValuesDotDotOld() const
240 : {
241 0 : return _element_data->dofValuesDotDotOld();
242 : }
243 :
244 : template <typename OutputType>
245 : const typename MooseVariableFV<OutputType>::DofValues &
246 0 : MooseVariableFV<OutputType>::dofValuesDotNeighbor() const
247 : {
248 0 : return _neighbor_data->dofValuesDot();
249 : }
250 :
251 : template <typename OutputType>
252 : const typename MooseVariableFV<OutputType>::DofValues &
253 0 : MooseVariableFV<OutputType>::dofValuesDotDotNeighbor() const
254 : {
255 0 : return _neighbor_data->dofValuesDotDot();
256 : }
257 :
258 : template <typename OutputType>
259 : const typename MooseVariableFV<OutputType>::DofValues &
260 0 : MooseVariableFV<OutputType>::dofValuesDotOldNeighbor() const
261 : {
262 0 : return _neighbor_data->dofValuesDotOld();
263 : }
264 :
265 : template <typename OutputType>
266 : const typename MooseVariableFV<OutputType>::DofValues &
267 0 : MooseVariableFV<OutputType>::dofValuesDotDotOldNeighbor() const
268 : {
269 0 : return _neighbor_data->dofValuesDotDotOld();
270 : }
271 :
272 : template <typename OutputType>
273 : const MooseArray<Number> &
274 0 : MooseVariableFV<OutputType>::dofValuesDuDotDu() const
275 : {
276 0 : return _element_data->dofValuesDuDotDu();
277 : }
278 :
279 : template <typename OutputType>
280 : const MooseArray<Number> &
281 0 : MooseVariableFV<OutputType>::dofValuesDuDotDotDu() const
282 : {
283 0 : return _element_data->dofValuesDuDotDotDu();
284 : }
285 :
286 : template <typename OutputType>
287 : const MooseArray<Number> &
288 0 : MooseVariableFV<OutputType>::dofValuesDuDotDuNeighbor() const
289 : {
290 0 : return _neighbor_data->dofValuesDuDotDu();
291 : }
292 :
293 : template <typename OutputType>
294 : const MooseArray<Number> &
295 0 : MooseVariableFV<OutputType>::dofValuesDuDotDotDuNeighbor() const
296 : {
297 0 : return _neighbor_data->dofValuesDuDotDotDu();
298 : }
299 :
300 : template <typename OutputType>
301 : void
302 201882 : MooseVariableFV<OutputType>::prepareIC()
303 : {
304 201882 : _element_data->prepareIC();
305 201882 : }
306 :
307 : template <typename OutputType>
308 : void
309 16235168 : MooseVariableFV<OutputType>::computeElemValues()
310 : {
311 16235168 : _element_data->setGeometry(Moose::Volume);
312 16235168 : _element_data->computeValues();
313 16235168 : }
314 :
315 : template <typename OutputType>
316 : void
317 190214 : MooseVariableFV<OutputType>::computeElemValuesFace()
318 : {
319 190214 : _element_data->setGeometry(Moose::Face);
320 190214 : _element_data->computeValues();
321 190214 : }
322 :
323 : template <typename OutputType>
324 : void
325 147180 : MooseVariableFV<OutputType>::computeNeighborValuesFace()
326 : {
327 147180 : _neighbor_data->setGeometry(Moose::Face);
328 147180 : _neighbor_data->computeValues();
329 147180 : }
330 :
331 : template <typename OutputType>
332 : void
333 0 : MooseVariableFV<OutputType>::computeNeighborValues()
334 : {
335 0 : _neighbor_data->setGeometry(Moose::Volume);
336 0 : _neighbor_data->computeValues();
337 0 : }
338 :
339 : template <typename OutputType>
340 : void
341 16151036 : MooseVariableFV<OutputType>::computeFaceValues(const FaceInfo & fi)
342 : {
343 16151036 : _element_data->setGeometry(Moose::Face);
344 16151036 : _neighbor_data->setGeometry(Moose::Face);
345 :
346 16151036 : const auto facetype = fi.faceType(std::make_pair(this->number(), this->sys().number()));
347 16151036 : if (facetype == FaceInfo::VarFaceNeighbors::NEITHER)
348 0 : return;
349 16151036 : else if (facetype == FaceInfo::VarFaceNeighbors::BOTH)
350 : {
351 14446481 : _element_data->computeValuesFace(fi);
352 14446481 : _neighbor_data->computeValuesFace(fi);
353 : }
354 1704555 : else if (facetype == FaceInfo::VarFaceNeighbors::ELEM)
355 1677842 : _element_data->computeValuesFace(fi);
356 26713 : else if (facetype == FaceInfo::VarFaceNeighbors::NEIGHBOR)
357 26713 : _neighbor_data->computeValuesFace(fi);
358 : else
359 0 : mooseError("robert wrote broken MooseVariableFV code");
360 : }
361 :
362 : template <typename OutputType>
363 : OutputType
364 0 : MooseVariableFV<OutputType>::getValue(const Elem * elem) const
365 : {
366 0 : Moose::initDofIndices(const_cast<MooseVariableFV<OutputType> &>(*this), *elem);
367 : mooseAssert(this->_dof_indices.size() == 1, "Wrong size for dof indices");
368 0 : OutputType value = (*this->_sys.currentSolution())(this->_dof_indices[0]);
369 0 : return value;
370 : }
371 :
372 : template <typename OutputType>
373 : typename OutputTools<OutputType>::OutputGradient
374 0 : MooseVariableFV<OutputType>::getGradient(const Elem * /*elem*/) const
375 : {
376 0 : return {};
377 : }
378 :
379 : template <typename OutputType>
380 : void
381 0 : MooseVariableFV<OutputType>::setNodalValue(const OutputType & /*value*/, unsigned int /*idx*/)
382 : {
383 0 : mooseError("FV variables do not support setNodalValue");
384 : }
385 :
386 : template <typename OutputType>
387 : void
388 278372 : MooseVariableFV<OutputType>::setDofValue(const DofValue & value, unsigned int index)
389 : {
390 278372 : _element_data->setDofValue(value, index);
391 278372 : }
392 :
393 : template <typename OutputType>
394 : void
395 4100 : MooseVariableFV<OutputType>::setDofValues(const DenseVector<DofValue> & values)
396 : {
397 4100 : _element_data->setDofValues(values);
398 4100 : }
399 :
400 : template <typename OutputType>
401 : void
402 0 : MooseVariableFV<OutputType>::setLowerDofValues(const DenseVector<DofValue> &)
403 : {
404 0 : lowerDError();
405 : }
406 :
407 : template <typename OutputType>
408 : std::pair<bool, const FVDirichletBCBase *>
409 217000507 : MooseVariableFV<OutputType>::getDirichletBC(const FaceInfo & fi) const
410 : {
411 220881282 : for (const auto bnd_id : fi.boundaryIDs())
412 14017311 : if (auto it = _boundary_id_to_dirichlet_bc.find(bnd_id);
413 14017311 : it != _boundary_id_to_dirichlet_bc.end())
414 10136536 : return {true, it->second};
415 :
416 206863971 : return {false, nullptr};
417 : }
418 :
419 : template <typename OutputType>
420 : std::pair<bool, std::vector<const FVFluxBC *>>
421 1737710 : MooseVariableFV<OutputType>::getFluxBCs(const FaceInfo & fi) const
422 : {
423 3064197 : for (const auto bnd_id : fi.boundaryIDs())
424 1862566 : if (auto it = _boundary_id_to_flux_bc.find(bnd_id); it != _boundary_id_to_flux_bc.end())
425 536079 : return {true, it->second};
426 :
427 1201631 : return std::make_pair(false, std::vector<const FVFluxBC *>());
428 : }
429 :
430 : template <typename OutputType>
431 : ADReal
432 206981309 : MooseVariableFV<OutputType>::getElemValue(const Elem * const elem, const StateArg & state) const
433 : {
434 : mooseAssert(elem,
435 : "The elem shall exist! This typically occurs when the "
436 : "user wants to evaluate non-existing elements (nullptr) at physical boundaries.");
437 : mooseAssert(
438 : this->hasBlocks(elem->subdomain_id()),
439 : "The variable should be defined on the element's subdomain! This typically occurs when the "
440 : "user wants to evaluate the elements right next to the boundary of two variables (block "
441 : "boundary). The subdomain which is queried: " +
442 : Moose::stringify(this->activeSubdomains()) + " the subdomain of the element " +
443 : std::to_string(elem->subdomain_id()));
444 :
445 206981309 : Moose::initDofIndices(const_cast<MooseVariableFV<OutputType> &>(*this), *elem);
446 :
447 : mooseAssert(
448 : this->_dof_indices.size() == 1,
449 : "There should only be one dof-index for a constant monomial variable on any given element");
450 :
451 206981309 : const dof_id_type index = this->_dof_indices[0];
452 :
453 : // It's not safe to use solutionState(0) because it returns the libMesh System solution member
454 : // which is wrong during things like finite difference Jacobian evaluation, e.g. when PETSc
455 : // perturbs the solution vector we feed these perturbations into the current_local_solution
456 : // while the libMesh solution is frozen in the non-perturbed state
457 206981309 : const auto & global_soln =
458 206981309 : (state.state == 0)
459 206981309 : ? *this->_sys.currentSolution()
460 145705940 : : std::as_const(this->_sys).solutionState(state.state, state.iteration_type);
461 :
462 206981309 : ADReal value = global_soln(index);
463 :
464 221980556 : if (ADReal::do_derivatives && state.state == 0 &&
465 14999247 : this->_sys.number() == this->_subproblem.currentNlSysNum())
466 14709228 : Moose::derivInsert(value.derivatives(), index, 1.);
467 :
468 206981309 : return value;
469 0 : }
470 :
471 : template <typename OutputType>
472 : bool
473 212408473 : MooseVariableFV<OutputType>::isDirichletBoundaryFace(const FaceInfo & fi,
474 : const Elem *,
475 : const Moose::StateArg &) const
476 : {
477 212408473 : const auto & pr = getDirichletBC(fi);
478 :
479 : // First member of this pair indicates whether we have a DirichletBC
480 212408473 : return pr.first;
481 : }
482 :
483 : template <typename OutputType>
484 : ADReal
485 3390403 : MooseVariableFV<OutputType>::getDirichletBoundaryFaceValue(const FaceInfo & fi,
486 : const Elem * const libmesh_dbg_var(elem),
487 : const Moose::StateArg & state) const
488 : {
489 : mooseAssert(isDirichletBoundaryFace(fi, elem, state),
490 : "This function should only be called on Dirichlet boundary faces.");
491 :
492 3390403 : const auto & diri_pr = getDirichletBC(fi);
493 :
494 : mooseAssert(diri_pr.first,
495 : "This functor should only be called if we are on a Dirichlet boundary face.");
496 :
497 3390403 : const FVDirichletBCBase & bc = *diri_pr.second;
498 :
499 6780806 : return ADReal(bc.boundaryValue(fi, state));
500 : }
501 :
502 : template <typename OutputType>
503 : bool
504 139556218 : MooseVariableFV<OutputType>::isExtrapolatedBoundaryFace(const FaceInfo & fi,
505 : const Elem * const elem,
506 : const Moose::StateArg & state) const
507 : {
508 139556218 : if (isDirichletBoundaryFace(fi, elem, state))
509 2519563 : return false;
510 : else
511 137036655 : return !this->isInternalFace(fi);
512 : }
513 :
514 : template <typename OutputType>
515 : ADReal
516 515893 : MooseVariableFV<OutputType>::getExtrapolatedBoundaryFaceValue(const FaceInfo & fi,
517 : const bool two_term_expansion,
518 : const bool correct_skewness,
519 : const Elem * elem_to_extrapolate_from,
520 : const StateArg & state) const
521 : {
522 : mooseAssert(
523 : isExtrapolatedBoundaryFace(fi, elem_to_extrapolate_from, state) || !two_term_expansion,
524 : "We allow Dirichlet boundary conditions to call this method. However, the only way to "
525 : "ensure we don't have infinite recursion, with Green Gauss gradients calling back to the "
526 : "Dirichlet boundary condition calling back to this method, is to do a one term expansion");
527 :
528 515893 : ADReal boundary_value;
529 : bool elem_to_extrapolate_from_is_fi_elem;
530 515893 : std::tie(elem_to_extrapolate_from, elem_to_extrapolate_from_is_fi_elem) =
531 0 : [this, &fi, elem_to_extrapolate_from]() -> std::pair<const Elem *, bool>
532 : {
533 515893 : if (elem_to_extrapolate_from)
534 : // Somebody already specified the element to extropolate from
535 510593 : return {elem_to_extrapolate_from, elem_to_extrapolate_from == fi.elemPtr()};
536 : else
537 : {
538 5300 : const auto [elem_guaranteed_to_have_dofs,
539 5300 : other_elem,
540 5300 : elem_guaranteed_to_have_dofs_is_fi_elem] =
541 : Moose::FV::determineElemOneAndTwo(fi, *this);
542 : // We only care about the element guaranteed to have degrees of freedom and current C++
543 : // doesn't allow us to not assign one of the returned items like python does
544 5300 : libmesh_ignore(other_elem);
545 : // We will extrapolate from the element guaranteed to have degrees of freedom
546 5300 : return {elem_guaranteed_to_have_dofs, elem_guaranteed_to_have_dofs_is_fi_elem};
547 : }
548 515893 : }();
549 :
550 515893 : if (two_term_expansion)
551 : {
552 30303 : const Point vector_to_face = elem_to_extrapolate_from_is_fi_elem
553 30303 : ? (fi.faceCentroid() - fi.elemCentroid())
554 1835 : : (fi.faceCentroid() - fi.neighborCentroid());
555 30303 : boundary_value = adGradSln(elem_to_extrapolate_from, state, correct_skewness) * vector_to_face +
556 : getElemValue(elem_to_extrapolate_from, state);
557 : }
558 : else
559 485590 : boundary_value = getElemValue(elem_to_extrapolate_from, state);
560 :
561 1031786 : return boundary_value;
562 0 : }
563 :
564 : template <typename OutputType>
565 : ADReal
566 359653 : MooseVariableFV<OutputType>::getBoundaryFaceValue(const FaceInfo & fi,
567 : const StateArg & state,
568 : const bool correct_skewness) const
569 : {
570 : mooseAssert(!this->isInternalFace(fi),
571 : "A boundary face value has been requested on an internal face.");
572 :
573 359653 : if (isDirichletBoundaryFace(fi, nullptr, state))
574 357151 : return getDirichletBoundaryFaceValue(fi, nullptr, state);
575 2502 : else if (isExtrapolatedBoundaryFace(fi, nullptr, state))
576 2502 : return getExtrapolatedBoundaryFaceValue(
577 2502 : fi, _two_term_boundary_expansion, correct_skewness, nullptr, state);
578 :
579 0 : mooseError("Unknown boundary face type!");
580 : }
581 :
582 : template <typename OutputType>
583 : const VectorValue<ADReal> &
584 29050005 : MooseVariableFV<OutputType>::adGradSln(const Elem * const elem,
585 : const StateArg & state,
586 : const bool correct_skewness) const
587 : {
588 : // We ensure that no caching takes place when we compute skewness-corrected
589 : // quantities.
590 29050005 : if (_cache_cell_gradients && !correct_skewness && state.state == 0)
591 : {
592 14511360 : auto it = _elem_to_grad.find(elem);
593 :
594 14511360 : if (it != _elem_to_grad.end())
595 10849186 : return it->second;
596 : }
597 :
598 18200819 : auto grad = FV::greenGaussGradient(
599 18200819 : ElemArg({elem, correct_skewness}), state, *this, _two_term_boundary_expansion, this->_mesh);
600 :
601 18200819 : if (_cache_cell_gradients && !correct_skewness && state.state == 0)
602 : {
603 3662174 : auto pr = _elem_to_grad.emplace(elem, std::move(grad));
604 : mooseAssert(pr.second, "Insertion should have just happened.");
605 3662174 : return pr.first->second;
606 : }
607 : else
608 : {
609 14538645 : _temp_cell_gradient = std::move(grad);
610 14538645 : return _temp_cell_gradient;
611 : }
612 18200819 : }
613 :
614 : template <typename OutputType>
615 : VectorValue<ADReal>
616 10153921 : MooseVariableFV<OutputType>::uncorrectedAdGradSln(const FaceInfo & fi,
617 : const StateArg & state,
618 : const bool correct_skewness) const
619 : {
620 10153921 : const auto face_type = fi.faceType(std::make_pair(this->number(), this->sys().number()));
621 : mooseAssert(face_type != FaceInfo::VarFaceNeighbors::NEITHER,
622 : "Gradient requested on a face where the variable is defined on neither side.");
623 :
624 10153921 : const bool var_defined_on_elem = (face_type == FaceInfo::VarFaceNeighbors::BOTH) ||
625 : (face_type == FaceInfo::VarFaceNeighbors::ELEM);
626 10153921 : const Elem * const elem_one = var_defined_on_elem ? &fi.elem() : fi.neighborPtr();
627 10153921 : const Elem * const elem_two = var_defined_on_elem ? fi.neighborPtr() : &fi.elem();
628 :
629 10153921 : const VectorValue<ADReal> elem_one_grad = adGradSln(elem_one, state, correct_skewness);
630 :
631 : // If we have a neighbor then we interpolate between the two to the face. If we do not, then we
632 : // apply a zero Hessian assumption and use the element centroid gradient as the uncorrected face
633 : // gradient
634 10153921 : if (face_type == FaceInfo::VarFaceNeighbors::BOTH)
635 : {
636 : mooseAssert(elem_two, "Face type indicates BOTH but neighbor information is missing.");
637 9843544 : const VectorValue<ADReal> & elem_two_grad = adGradSln(elem_two, state, correct_skewness);
638 :
639 : // Uncorrected gradient value
640 9843544 : return Moose::FV::linearInterpolation(elem_one_grad, elem_two_grad, fi, var_defined_on_elem);
641 : }
642 : else
643 310377 : return elem_one_grad;
644 10153921 : }
645 :
646 : template <typename OutputType>
647 : VectorValue<ADReal>
648 13646794 : MooseVariableFV<OutputType>::adGradSln(const FaceInfo & fi,
649 : const StateArg & state,
650 : const bool correct_skewness) const
651 : {
652 13646794 : const bool var_defined_on_elem = this->hasBlocks(fi.elem().subdomain_id());
653 13646794 : const Elem * const elem = &fi.elem();
654 13646794 : const Elem * const neighbor = fi.neighborPtr();
655 :
656 13646794 : const bool is_internal_face = this->isInternalFace(fi);
657 :
658 13646794 : const ADReal side_one_value = (!is_internal_face && !var_defined_on_elem)
659 13646794 : ? getBoundaryFaceValue(fi, state, correct_skewness)
660 : : getElemValue(elem, state);
661 13646794 : const ADReal side_two_value = (var_defined_on_elem && !is_internal_face)
662 13646794 : ? getBoundaryFaceValue(fi, state, correct_skewness)
663 : : getElemValue(neighbor, state);
664 :
665 13646794 : const auto delta =
666 13646794 : this->isInternalFace(fi)
667 13646794 : ? fi.dCNMag()
668 358373 : : (fi.faceCentroid() - (var_defined_on_elem ? fi.elemCentroid() : fi.neighborCentroid()))
669 358373 : .norm();
670 :
671 : // This is the component of the gradient which is parallel to the line connecting
672 : // the cell centers. Therefore, we can use our second order, central difference
673 : // scheme to approximate it.
674 13646794 : auto face_grad = ((side_two_value - side_one_value) / delta) * fi.eCN();
675 :
676 : // We only need non-orthogonal correctors in 2+ dimensions
677 13646794 : if (this->_mesh.dimension() > 1)
678 : {
679 : // We are using an orthogonal approach for the non-orthogonal correction, for more information
680 : // see Hrvoje Jasak's PhD Thesis (Imperial College, 1996)
681 10153921 : const auto & interpolated_gradient = uncorrectedAdGradSln(fi, state, correct_skewness);
682 10153921 : face_grad += interpolated_gradient - (interpolated_gradient * fi.eCN()) * fi.eCN();
683 10153921 : }
684 :
685 27293588 : return face_grad;
686 13646794 : }
687 :
688 : template <typename OutputType>
689 : void
690 111546 : MooseVariableFV<OutputType>::residualSetup()
691 : {
692 111546 : if (!_dirichlet_map_setup)
693 16438 : determineBoundaryToDirichletBCMap();
694 111546 : if (!_flux_map_setup)
695 16438 : determineBoundaryToFluxBCMap();
696 :
697 111546 : clearCaches();
698 111546 : }
699 :
700 : template <typename OutputType>
701 : void
702 35711 : MooseVariableFV<OutputType>::jacobianSetup()
703 : {
704 35711 : clearCaches();
705 35711 : }
706 :
707 : template <typename OutputType>
708 : void
709 147257 : MooseVariableFV<OutputType>::clearCaches()
710 : {
711 147257 : _elem_to_grad.clear();
712 147257 : }
713 :
714 : template <typename OutputType>
715 : unsigned int
716 4893 : MooseVariableFV<OutputType>::oldestSolutionStateRequested() const
717 : {
718 4893 : unsigned int state = 0;
719 4893 : state = std::max(state, _element_data->oldestSolutionStateRequested());
720 4893 : state = std::max(state, _neighbor_data->oldestSolutionStateRequested());
721 4893 : return state;
722 : }
723 :
724 : template <typename OutputType>
725 : void
726 27510 : MooseVariableFV<OutputType>::clearAllDofIndices()
727 : {
728 27510 : _element_data->clearDofIndices();
729 27510 : _neighbor_data->clearDofIndices();
730 27510 : }
731 :
732 : template <typename OutputType>
733 : typename MooseVariableFV<OutputType>::ValueType
734 72492602 : MooseVariableFV<OutputType>::evaluate(const FaceArg & face, const StateArg & state) const
735 : {
736 72492602 : const FaceInfo * const fi = face.fi;
737 : mooseAssert(fi, "The face information must be non-null");
738 72492602 : if (isDirichletBoundaryFace(*fi, face.face_side, state))
739 3033252 : return getDirichletBoundaryFaceValue(*fi, face.face_side, state);
740 69459350 : else if (isExtrapolatedBoundaryFace(*fi, face.face_side, state))
741 : {
742 513391 : bool two_term_boundary_expansion = _two_term_boundary_expansion;
743 513391 : if (face.limiter_type == Moose::FV::LimiterType::Upwind)
744 135 : if ((face.elem_is_upwind && face.face_side == fi->elemPtr()) ||
745 0 : (!face.elem_is_upwind && face.face_side == fi->neighborPtr()))
746 135 : two_term_boundary_expansion = false;
747 513391 : return getExtrapolatedBoundaryFaceValue(
748 513391 : *fi, two_term_boundary_expansion, face.correct_skewness, face.face_side, state);
749 : }
750 : else
751 : {
752 : mooseAssert(this->isInternalFace(*fi),
753 : "We must be either Dirichlet, extrapolated, or internal");
754 68945959 : return Moose::FV::interpolate(*this, face, state);
755 : }
756 : }
757 :
758 : template <typename OutputType>
759 : typename MooseVariableFV<OutputType>::ValueType
760 0 : MooseVariableFV<OutputType>::evaluate(const NodeArg & node_arg, const StateArg & state) const
761 : {
762 0 : const auto & node_to_elem_map = this->_mesh.nodeToElemMap();
763 0 : const auto & elem_ids = libmesh_map_find(node_to_elem_map, node_arg.node->id());
764 0 : ValueType sum = 0;
765 0 : Real total_weight = 0;
766 : mooseAssert(elem_ids.size(), "There should always be at least one element connected to a node");
767 0 : for (const auto elem_id : elem_ids)
768 : {
769 0 : const Elem * const elem = this->_mesh.queryElemPtr(elem_id);
770 : mooseAssert(elem, "We should have this element available");
771 0 : if (!this->hasBlocks(elem->subdomain_id()))
772 0 : continue;
773 0 : const ElemPointArg elem_point{
774 0 : elem, *node_arg.node, _face_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage};
775 0 : const auto weight = 1 / (*node_arg.node - elem->vertex_average()).norm();
776 0 : sum += weight * (*this)(elem_point, state);
777 0 : total_weight += weight;
778 : }
779 0 : return sum / total_weight;
780 0 : }
781 :
782 : template <typename OutputType>
783 : typename MooseVariableFV<OutputType>::DotType
784 : MooseVariableFV<OutputType>::evaluateDot(const ElemArg &, const StateArg &) const
785 : {
786 : mooseError("evaluateDot not implemented for this class of finite volume variables");
787 : }
788 :
789 : template <>
790 : ADReal
791 56292 : MooseVariableFV<Real>::evaluateDot(const ElemArg & elem_arg, const StateArg & state) const
792 : {
793 56292 : const Elem * const elem = elem_arg.elem;
794 : mooseAssert(state.state == 0,
795 : "We dot not currently support any time derivative evaluations other than for the "
796 : "current time-step");
797 : mooseAssert(_time_integrator && _time_integrator->dt(),
798 : "A time derivative is being requested but we do not have a time integrator so we'll "
799 : "have no idea how to compute it");
800 :
801 56292 : Moose::initDofIndices(const_cast<MooseVariableFV<Real> &>(*this), *elem);
802 :
803 : mooseAssert(
804 : this->_dof_indices.size() == 1,
805 : "There should only be one dof-index for a constant monomial variable on any given element");
806 :
807 56292 : const dof_id_type dof_index = this->_dof_indices[0];
808 :
809 56292 : if (_var_kind == Moose::VAR_SOLVER)
810 : {
811 56292 : ADReal dot = (*_solution)(dof_index);
812 105204 : if (ADReal::do_derivatives && state.state == 0 &&
813 48912 : _sys.number() == _subproblem.currentNlSysNum())
814 48912 : Moose::derivInsert(dot.derivatives(), dof_index, 1.);
815 56292 : _time_integrator->computeADTimeDerivatives(dot, dof_index, _ad_real_dummy);
816 56292 : return dot;
817 56292 : }
818 : else
819 0 : return (*_sys.solutionUDot())(dof_index);
820 : }
821 :
822 : template <>
823 : ADReal
824 0 : MooseVariableFV<Real>::evaluateDot(const FaceArg & face, const StateArg & state) const
825 : {
826 0 : const FaceInfo * const fi = face.fi;
827 : mooseAssert(fi, "The face information must be non-null");
828 0 : if (isDirichletBoundaryFace(*fi, face.face_side, state))
829 0 : return ADReal(0.0); // No time derivative if boundary value is set
830 0 : else if (isExtrapolatedBoundaryFace(*fi, face.face_side, state))
831 : {
832 : mooseAssert(face.face_side && this->hasBlocks(face.face_side->subdomain_id()),
833 : "If we are an extrapolated boundary face, then our FunctorBase::checkFace method "
834 : "should have assigned a non-null element that we are defined on");
835 0 : const auto elem_arg = ElemArg({face.face_side, face.correct_skewness});
836 : // For extrapolated boundary faces, note that we take the value of the time derivative at the
837 : // cell in contact with the face
838 0 : return evaluateDot(elem_arg, state);
839 : }
840 : else
841 : {
842 : mooseAssert(this->isInternalFace(*fi),
843 : "We must be either Dirichlet, extrapolated, or internal");
844 0 : return Moose::FV::interpolate<ADReal, FunctorEvaluationKind::Dot>(*this, face, state);
845 : }
846 : }
847 :
848 : template <>
849 : ADReal
850 0 : MooseVariableFV<Real>::evaluateDot(const ElemQpArg & elem_qp, const StateArg & state) const
851 : {
852 0 : return evaluateDot(ElemArg({elem_qp.elem, /*correct_skewness*/ false}), state);
853 : }
854 :
855 : template <typename OutputType>
856 : void
857 108528 : MooseVariableFV<OutputType>::prepareAux()
858 : {
859 108528 : _element_data->prepareAux();
860 108528 : _neighbor_data->prepareAux();
861 108528 : }
862 :
863 : template <typename OutputType>
864 : void
865 21683 : MooseVariableFV<OutputType>::determineBoundaryToDirichletBCMap()
866 : {
867 : mooseAssert(!Threads::in_threads,
868 : "This routine has not been implemented for threads. Please query this routine before "
869 : "a threaded region or contact a MOOSE developer to discuss.");
870 :
871 21683 : _boundary_id_to_dirichlet_bc.clear();
872 21683 : std::vector<FVDirichletBCBase *> bcs;
873 :
874 : // I believe because query() returns by value but condition returns by reference that binding to a
875 : // const lvalue reference results in the query() getting destructed and us holding onto a dangling
876 : // reference. I think that condition returned by value we would be able to bind to a const lvalue
877 : // reference here. But as it is we'll bind to a regular lvalue
878 21683 : const auto base_query = this->_subproblem.getMooseApp()
879 21683 : .theWarehouse()
880 : .query()
881 21683 : .template condition<AttribSystem>("FVDirichletBC")
882 21683 : .template condition<AttribThread>(_tid)
883 21683 : .template condition<AttribVar>(_var_num)
884 21683 : .template condition<AttribSysNum>(this->_sys.number());
885 :
886 102650 : for (const auto bnd_id : this->_mesh.getBoundaryIDs())
887 : {
888 80967 : auto base_query_copy = base_query;
889 161934 : base_query_copy.template condition<AttribBoundaries>(std::set<BoundaryID>({bnd_id}))
890 80967 : .queryInto(bcs);
891 : mooseAssert(bcs.size() <= 1, "cannot have multiple dirichlet BCs on the same boundary");
892 80967 : if (!bcs.empty())
893 32226 : _boundary_id_to_dirichlet_bc.emplace(bnd_id, bcs[0]);
894 : }
895 :
896 21683 : _dirichlet_map_setup = true;
897 21683 : }
898 :
899 : template <typename OutputType>
900 : void
901 21683 : MooseVariableFV<OutputType>::determineBoundaryToFluxBCMap()
902 : {
903 : mooseAssert(!Threads::in_threads,
904 : "This routine has not been implemented for threads. Please query this routine before "
905 : "a threaded region or contact a MOOSE developer to discuss.");
906 :
907 21683 : _boundary_id_to_flux_bc.clear();
908 21683 : std::vector<const FVFluxBC *> bcs;
909 :
910 : // I believe because query() returns by value but condition returns by reference that binding to a
911 : // const lvalue reference results in the query() getting destructed and us holding onto a dangling
912 : // reference. I think that condition returned by value we would be able to bind to a const lvalue
913 : // reference here. But as it is we'll bind to a regular lvalue
914 21683 : const auto base_query = this->_subproblem.getMooseApp()
915 21683 : .theWarehouse()
916 : .query()
917 21683 : .template condition<AttribSystem>("FVFluxBC")
918 21683 : .template condition<AttribThread>(_tid)
919 21683 : .template condition<AttribVar>(_var_num)
920 21683 : .template condition<AttribSysNum>(this->_sys.number());
921 :
922 102650 : for (const auto bnd_id : this->_mesh.getBoundaryIDs())
923 : {
924 80967 : auto base_query_copy = base_query;
925 161934 : base_query_copy.template condition<AttribBoundaries>(std::set<BoundaryID>({bnd_id}))
926 80967 : .queryInto(bcs);
927 80967 : if (!bcs.empty())
928 21896 : _boundary_id_to_flux_bc.emplace(bnd_id, bcs);
929 : }
930 :
931 21683 : _flux_map_setup = true;
932 21683 : }
933 :
934 : template <typename OutputType>
935 : void
936 5391 : MooseVariableFV<OutputType>::sizeMatrixTagData()
937 : {
938 5391 : _element_data->sizeMatrixTagData();
939 5391 : _neighbor_data->sizeMatrixTagData();
940 5391 : }
941 :
942 : template class MooseVariableFV<Real>;
943 : // TODO: implement vector fv variable support. This will require some template
944 : // specializations for various member functions in this and the FV variable
945 : // classes. And then you will need to uncomment out the line below:
946 : // template class MooseVariableFV<RealVectorValue>;
|
