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 25064 : MooseVariableFV<OutputType>::validParams()
35 : {
36 25064 : InputParameters params = MooseVariableField<OutputType>::validParams();
37 25064 : params.set<bool>("fv") = true;
38 25064 : params.set<MooseEnum>("family") = "MONOMIAL";
39 25064 : params.set<MooseEnum>("order") = "CONSTANT";
40 75192 : params.template addParam<bool>(
41 : "two_term_boundary_expansion",
42 50128 : 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 25064 : MooseEnum face_interp_method("average skewness-corrected", "average");
49 25064 : params.template addParam<MooseEnum>("face_interp_method",
50 : face_interp_method,
51 : "Switch that can select between face interpolation methods.");
52 75192 : params.template addParam<bool>(
53 50128 : "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 25064 : 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 12731 : unsigned short layers = 1;
63 12731 : if (obj_params.get<MooseEnum>("face_interp_method") == "skewness-corrected")
64 420 : layers = 2;
65 :
66 12731 : rm_params.set<unsigned short>("layers") = layers;
67 : });
68 50128 : return params;
69 25064 : }
70 :
71 : template <typename OutputType>
72 6374 : MooseVariableFV<OutputType>::MooseVariableFV(const InputParameters & parameters)
73 : : MooseVariableField<OutputType>(parameters),
74 6374 : _solution(this->_sys.currentSolution()),
75 6374 : _phi(this->_assembly.template fePhi<OutputShape>(FEType(CONSTANT, MONOMIAL))),
76 6374 : _grad_phi(this->_assembly.template feGradPhi<OutputShape>(FEType(CONSTANT, MONOMIAL))),
77 6374 : _phi_face(this->_assembly.template fePhiFace<OutputShape>(FEType(CONSTANT, MONOMIAL))),
78 6374 : _grad_phi_face(this->_assembly.template feGradPhiFace<OutputShape>(FEType(CONSTANT, MONOMIAL))),
79 6374 : _phi_face_neighbor(
80 6374 : this->_assembly.template fePhiFaceNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
81 6374 : _grad_phi_face_neighbor(
82 6374 : this->_assembly.template feGradPhiFaceNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
83 6374 : _phi_neighbor(this->_assembly.template fePhiNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
84 6374 : _grad_phi_neighbor(
85 6374 : this->_assembly.template feGradPhiNeighbor<OutputShape>(FEType(CONSTANT, MONOMIAL))),
86 6374 : _prev_elem(nullptr),
87 12748 : _two_term_boundary_expansion(this->isParamValid("two_term_boundary_expansion")
88 11074 : ? this->template getParam<bool>("two_term_boundary_expansion")
89 : : true),
90 12748 : _cache_cell_gradients(this->isParamValid("cache_cell_gradients")
91 11074 : ? this->template getParam<bool>("cache_cell_gradients")
92 19122 : : true)
93 : {
94 6374 : _element_data = std::make_unique<MooseVariableDataFV<OutputType>>(
95 6374 : *this, _sys, _tid, Moose::ElementType::Element, this->_assembly.elem());
96 6374 : _neighbor_data = std::make_unique<MooseVariableDataFV<OutputType>>(
97 6374 : *this, _sys, _tid, Moose::ElementType::Neighbor, this->_assembly.neighbor());
98 :
99 6374 : if (this->isParamValid("face_interp_method"))
100 : {
101 4700 : const auto & interp_method = this->template getParam<MooseEnum>("face_interp_method");
102 4700 : if (interp_method == "average")
103 4560 : _face_interp_method = Moose::FV::InterpMethod::Average;
104 140 : else if (interp_method == "skewness-corrected")
105 140 : _face_interp_method = Moose::FV::InterpMethod::SkewCorrectedAverage;
106 : }
107 : else
108 1674 : _face_interp_method = Moose::FV::InterpMethod::Average;
109 15774 : }
110 :
111 : template <typename OutputType>
112 : void
113 20966107 : MooseVariableFV<OutputType>::clearDofIndices()
114 : {
115 20966107 : _element_data->clearDofIndices();
116 20966107 : }
117 :
118 : template <typename OutputType>
119 : typename MooseVariableFV<OutputType>::OutputData
120 2326 : MooseVariableFV<OutputType>::getElementalValue(const Elem * elem, unsigned int idx) const
121 : {
122 2326 : return _element_data->getElementalValue(elem, Moose::Current, idx);
123 : }
124 :
125 : template <typename OutputType>
126 : typename MooseVariableFV<OutputType>::OutputData
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>::OutputData
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 430736 : MooseVariableFV<OutputType>::insert(NumericVector<Number> & residual)
142 : {
143 430736 : _element_data->insert(residual);
144 430736 : }
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>::DoFValue &
162 65 : MooseVariableFV<OutputType>::dofValues() const
163 : {
164 65 : return _element_data->dofValues();
165 : }
166 :
167 : template <typename OutputType>
168 : const typename MooseVariableFV<OutputType>::DoFValue &
169 0 : MooseVariableFV<OutputType>::dofValuesOld() const
170 : {
171 0 : return _element_data->dofValuesOld();
172 : }
173 :
174 : template <typename OutputType>
175 : const typename MooseVariableFV<OutputType>::DoFValue &
176 0 : MooseVariableFV<OutputType>::dofValuesOlder() const
177 : {
178 0 : return _element_data->dofValuesOlder();
179 : }
180 :
181 : template <typename OutputType>
182 : const typename MooseVariableFV<OutputType>::DoFValue &
183 0 : MooseVariableFV<OutputType>::dofValuesPreviousNL() const
184 : {
185 0 : return _element_data->dofValuesPreviousNL();
186 : }
187 :
188 : template <typename OutputType>
189 : const typename MooseVariableFV<OutputType>::DoFValue &
190 0 : MooseVariableFV<OutputType>::dofValuesNeighbor() const
191 : {
192 0 : return _neighbor_data->dofValues();
193 : }
194 :
195 : template <typename OutputType>
196 : const typename MooseVariableFV<OutputType>::DoFValue &
197 0 : MooseVariableFV<OutputType>::dofValuesOldNeighbor() const
198 : {
199 0 : return _neighbor_data->dofValuesOld();
200 : }
201 :
202 : template <typename OutputType>
203 : const typename MooseVariableFV<OutputType>::DoFValue &
204 0 : MooseVariableFV<OutputType>::dofValuesOlderNeighbor() const
205 : {
206 0 : return _neighbor_data->dofValuesOlder();
207 : }
208 :
209 : template <typename OutputType>
210 : const typename MooseVariableFV<OutputType>::DoFValue &
211 0 : MooseVariableFV<OutputType>::dofValuesPreviousNLNeighbor() const
212 : {
213 0 : return _neighbor_data->dofValuesPreviousNL();
214 : }
215 :
216 : template <typename OutputType>
217 : const typename MooseVariableFV<OutputType>::DoFValue &
218 0 : MooseVariableFV<OutputType>::dofValuesDot() const
219 : {
220 0 : return _element_data->dofValuesDot();
221 : }
222 :
223 : template <typename OutputType>
224 : const typename MooseVariableFV<OutputType>::DoFValue &
225 0 : MooseVariableFV<OutputType>::dofValuesDotDot() const
226 : {
227 0 : return _element_data->dofValuesDotDot();
228 : }
229 :
230 : template <typename OutputType>
231 : const typename MooseVariableFV<OutputType>::DoFValue &
232 0 : MooseVariableFV<OutputType>::dofValuesDotOld() const
233 : {
234 0 : return _element_data->dofValuesDotOld();
235 : }
236 :
237 : template <typename OutputType>
238 : const typename MooseVariableFV<OutputType>::DoFValue &
239 0 : MooseVariableFV<OutputType>::dofValuesDotDotOld() const
240 : {
241 0 : return _element_data->dofValuesDotDotOld();
242 : }
243 :
244 : template <typename OutputType>
245 : const typename MooseVariableFV<OutputType>::DoFValue &
246 0 : MooseVariableFV<OutputType>::dofValuesDotNeighbor() const
247 : {
248 0 : return _neighbor_data->dofValuesDot();
249 : }
250 :
251 : template <typename OutputType>
252 : const typename MooseVariableFV<OutputType>::DoFValue &
253 0 : MooseVariableFV<OutputType>::dofValuesDotDotNeighbor() const
254 : {
255 0 : return _neighbor_data->dofValuesDotDot();
256 : }
257 :
258 : template <typename OutputType>
259 : const typename MooseVariableFV<OutputType>::DoFValue &
260 0 : MooseVariableFV<OutputType>::dofValuesDotOldNeighbor() const
261 : {
262 0 : return _neighbor_data->dofValuesDotOld();
263 : }
264 :
265 : template <typename OutputType>
266 : const typename MooseVariableFV<OutputType>::DoFValue &
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 385530 : MooseVariableFV<OutputType>::prepareIC()
303 : {
304 385530 : _element_data->prepareIC();
305 385530 : }
306 :
307 : template <typename OutputType>
308 : void
309 26958530 : MooseVariableFV<OutputType>::computeElemValues()
310 : {
311 26958530 : _element_data->setGeometry(Moose::Volume);
312 26958530 : _element_data->computeValues();
313 26958530 : }
314 :
315 : template <typename OutputType>
316 : void
317 192068 : MooseVariableFV<OutputType>::computeElemValuesFace()
318 : {
319 192068 : _element_data->setGeometry(Moose::Face);
320 192068 : _element_data->computeValues();
321 192068 : }
322 :
323 : template <typename OutputType>
324 : void
325 152188 : MooseVariableFV<OutputType>::computeNeighborValuesFace()
326 : {
327 152188 : _neighbor_data->setGeometry(Moose::Face);
328 152188 : _neighbor_data->computeValues();
329 152188 : }
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 27086703 : MooseVariableFV<OutputType>::computeFaceValues(const FaceInfo & fi)
342 : {
343 27086703 : _element_data->setGeometry(Moose::Face);
344 27086703 : _neighbor_data->setGeometry(Moose::Face);
345 :
346 27086703 : const auto facetype = fi.faceType(std::make_pair(this->number(), this->sys().number()));
347 27086703 : if (facetype == FaceInfo::VarFaceNeighbors::NEITHER)
348 0 : return;
349 27086703 : else if (facetype == FaceInfo::VarFaceNeighbors::BOTH)
350 : {
351 24974221 : _element_data->computeValuesFace(fi);
352 24974221 : _neighbor_data->computeValuesFace(fi);
353 : }
354 2112482 : else if (facetype == FaceInfo::VarFaceNeighbors::ELEM)
355 2084964 : _element_data->computeValuesFace(fi);
356 27518 : else if (facetype == FaceInfo::VarFaceNeighbors::NEIGHBOR)
357 27518 : _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 390600 : MooseVariableFV<OutputType>::setDofValue(const OutputData & value, unsigned int index)
389 : {
390 390600 : _element_data->setDofValue(value, index);
391 390600 : }
392 :
393 : template <typename OutputType>
394 : void
395 4000 : MooseVariableFV<OutputType>::setDofValues(const DenseVector<OutputData> & values)
396 : {
397 4000 : _element_data->setDofValues(values);
398 4000 : }
399 :
400 : template <typename OutputType>
401 : void
402 0 : MooseVariableFV<OutputType>::setLowerDofValues(const DenseVector<OutputData> &)
403 : {
404 0 : lowerDError();
405 : }
406 :
407 : template <typename OutputType>
408 : std::pair<bool, const FVDirichletBCBase *>
409 252631446 : MooseVariableFV<OutputType>::getDirichletBC(const FaceInfo & fi) const
410 : {
411 257058338 : for (const auto bnd_id : fi.boundaryIDs())
412 15587941 : if (auto it = _boundary_id_to_dirichlet_bc.find(bnd_id);
413 15587941 : it != _boundary_id_to_dirichlet_bc.end())
414 11161049 : return {true, it->second};
415 :
416 241470397 : return {false, nullptr};
417 : }
418 :
419 : template <typename OutputType>
420 : std::pair<bool, std::vector<const FVFluxBC *>>
421 2165368 : MooseVariableFV<OutputType>::getFluxBCs(const FaceInfo & fi) const
422 : {
423 3767413 : for (const auto bnd_id : fi.boundaryIDs())
424 2312093 : if (auto it = _boundary_id_to_flux_bc.find(bnd_id); it != _boundary_id_to_flux_bc.end())
425 710048 : return {true, it->second};
426 :
427 2910640 : return std::make_pair(false, std::vector<const FVFluxBC *>());
428 : }
429 :
430 : template <typename OutputType>
431 : ADReal
432 243049329 : 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 243049329 : 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 243049329 : 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 243049329 : const auto & global_soln =
458 243049329 : (state.state == 0)
459 243049329 : ? *this->_sys.currentSolution()
460 145593984 : : std::as_const(this->_sys).solutionState(state.state, state.iteration_type);
461 :
462 243049329 : ADReal value = global_soln(index);
463 :
464 268513537 : if (ADReal::do_derivatives && state.state == 0 &&
465 25464208 : this->_sys.number() == this->_subproblem.currentNlSysNum())
466 25168687 : Moose::derivInsert(value.derivatives(), index, 1.);
467 :
468 243049329 : return value;
469 0 : }
470 :
471 : template <typename OutputType>
472 : bool
473 247443481 : MooseVariableFV<OutputType>::isDirichletBoundaryFace(const FaceInfo & fi,
474 : const Elem *,
475 : const Moose::StateArg &) const
476 : {
477 247443481 : const auto & pr = getDirichletBC(fi);
478 :
479 : // First member of this pair indicates whether we have a DirichletBC
480 247443481 : return pr.first;
481 : }
482 :
483 : template <typename OutputType>
484 : ADReal
485 3732645 : 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 3732645 : 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 3732645 : const FVDirichletBCBase & bc = *diri_pr.second;
498 :
499 7465290 : return ADReal(bc.boundaryValue(fi, state));
500 : }
501 :
502 : template <typename OutputType>
503 : bool
504 160774383 : MooseVariableFV<OutputType>::isExtrapolatedBoundaryFace(const FaceInfo & fi,
505 : const Elem * const elem,
506 : const Moose::StateArg & state) const
507 : {
508 160774383 : if (isDirichletBoundaryFace(fi, elem, state))
509 2611735 : return false;
510 : else
511 158162648 : return !this->isInternalFace(fi);
512 : }
513 :
514 : template <typename OutputType>
515 : ADReal
516 680730 : 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 680730 : ADReal boundary_value;
529 : bool elem_to_extrapolate_from_is_fi_elem;
530 680730 : 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 680730 : if (elem_to_extrapolate_from)
534 : // Somebody already specified the element to extropolate from
535 676612 : return {elem_to_extrapolate_from, elem_to_extrapolate_from == fi.elemPtr()};
536 : else
537 : {
538 4118 : const auto [elem_guaranteed_to_have_dofs,
539 4118 : other_elem,
540 4118 : 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 4118 : libmesh_ignore(other_elem);
545 : // We will extrapolate from the element guaranteed to have degrees of freedom
546 4118 : return {elem_guaranteed_to_have_dofs, elem_guaranteed_to_have_dofs_is_fi_elem};
547 : }
548 680730 : }();
549 :
550 680730 : if (two_term_expansion)
551 : {
552 114364 : const Point vector_to_face = elem_to_extrapolate_from_is_fi_elem
553 114364 : ? (fi.faceCentroid() - fi.elemCentroid())
554 1835 : : (fi.faceCentroid() - fi.neighborCentroid());
555 114364 : boundary_value = adGradSln(elem_to_extrapolate_from, state, correct_skewness) * vector_to_face +
556 : getElemValue(elem_to_extrapolate_from, state);
557 : }
558 : else
559 566366 : boundary_value = getElemValue(elem_to_extrapolate_from, state);
560 :
561 1361460 : return boundary_value;
562 0 : }
563 :
564 : template <typename OutputType>
565 : ADReal
566 428564 : 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 428564 : if (isDirichletBoundaryFace(fi, nullptr, state))
574 426062 : 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 38047849 : 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 38047849 : if (_cache_cell_gradients && !correct_skewness && state.state == 0)
591 : {
592 22636286 : auto it = _elem_to_grad.find(elem);
593 :
594 22636286 : if (it != _elem_to_grad.end())
595 17039239 : return it->second;
596 : }
597 :
598 21008610 : auto grad = FV::greenGaussGradient(
599 21008610 : ElemArg({elem, correct_skewness}), state, *this, _two_term_boundary_expansion, this->_mesh);
600 :
601 21008610 : if (_cache_cell_gradients && !correct_skewness && state.state == 0)
602 : {
603 5597047 : auto pr = _elem_to_grad.emplace(elem, std::move(grad));
604 : mooseAssert(pr.second, "Insertion should have just happened.");
605 5597047 : return pr.first->second;
606 : }
607 : else
608 : {
609 15411563 : _temp_cell_gradient = std::move(grad);
610 15411563 : return _temp_cell_gradient;
611 : }
612 21008610 : }
613 :
614 : template <typename OutputType>
615 : VectorValue<ADReal>
616 11877055 : MooseVariableFV<OutputType>::uncorrectedAdGradSln(const FaceInfo & fi,
617 : const StateArg & state,
618 : const bool correct_skewness) const
619 : {
620 11877055 : const bool var_defined_on_elem = this->hasBlocks(fi.elem().subdomain_id());
621 11877055 : const Elem * const elem_one = var_defined_on_elem ? &fi.elem() : fi.neighborPtr();
622 11877055 : const Elem * const elem_two = var_defined_on_elem ? fi.neighborPtr() : &fi.elem();
623 :
624 11877055 : const VectorValue<ADReal> elem_one_grad = adGradSln(elem_one, state, correct_skewness);
625 :
626 : // If we have a neighbor then we interpolate between the two to the face. If we do not, then we
627 : // apply a zero Hessian assumption and use the element centroid gradient as the uncorrected face
628 : // gradient
629 11877055 : if (elem_two && this->hasBlocks(elem_two->subdomain_id()))
630 : {
631 11502731 : const VectorValue<ADReal> & elem_two_grad = adGradSln(elem_two, state, correct_skewness);
632 :
633 : // Uncorrected gradient value
634 11502731 : return Moose::FV::linearInterpolation(elem_one_grad, elem_two_grad, fi, var_defined_on_elem);
635 : }
636 : else
637 374324 : return elem_one_grad;
638 11877055 : }
639 :
640 : template <typename OutputType>
641 : VectorValue<ADReal>
642 15508921 : MooseVariableFV<OutputType>::adGradSln(const FaceInfo & fi,
643 : const StateArg & state,
644 : const bool correct_skewness) const
645 : {
646 15508921 : const bool var_defined_on_elem = this->hasBlocks(fi.elem().subdomain_id());
647 15508921 : const Elem * const elem = &fi.elem();
648 15508921 : const Elem * const neighbor = fi.neighborPtr();
649 :
650 15508921 : const bool is_internal_face = this->isInternalFace(fi);
651 :
652 15508921 : const ADReal side_one_value = (!is_internal_face && !var_defined_on_elem)
653 15508921 : ? getBoundaryFaceValue(fi, state, correct_skewness)
654 : : getElemValue(elem, state);
655 15508921 : const ADReal side_two_value = (var_defined_on_elem && !is_internal_face)
656 15508921 : ? getBoundaryFaceValue(fi, state, correct_skewness)
657 : : getElemValue(neighbor, state);
658 :
659 15508921 : const auto delta =
660 15508921 : this->isInternalFace(fi)
661 15508921 : ? fi.dCNMag()
662 427284 : : (fi.faceCentroid() - (var_defined_on_elem ? fi.elemCentroid() : fi.neighborCentroid()))
663 427284 : .norm();
664 :
665 : // This is the component of the gradient which is parallel to the line connecting
666 : // the cell centers. Therefore, we can use our second order, central difference
667 : // scheme to approximate it.
668 15508921 : auto face_grad = ((side_two_value - side_one_value) / delta) * fi.eCN();
669 :
670 : // We only need non-orthogonal correctors in 2+ dimensions
671 15508921 : if (this->_mesh.dimension() > 1)
672 : {
673 : // We are using an orthogonal approach for the non-orthogonal correction, for more information
674 : // see Hrvoje Jasak's PhD Thesis (Imperial College, 1996)
675 11877055 : const auto & interpolated_gradient = uncorrectedAdGradSln(fi, state, correct_skewness);
676 11877055 : face_grad += interpolated_gradient - (interpolated_gradient * fi.eCN()) * fi.eCN();
677 11877055 : }
678 :
679 31017842 : return face_grad;
680 15508921 : }
681 :
682 : template <typename OutputType>
683 : void
684 259347 : MooseVariableFV<OutputType>::residualSetup()
685 : {
686 259347 : if (!_dirichlet_map_setup)
687 20881 : determineBoundaryToDirichletBCMap();
688 259347 : if (!_flux_map_setup)
689 20881 : determineBoundaryToFluxBCMap();
690 :
691 259347 : clearCaches();
692 259347 : }
693 :
694 : template <typename OutputType>
695 : void
696 57433 : MooseVariableFV<OutputType>::jacobianSetup()
697 : {
698 57433 : clearCaches();
699 57433 : }
700 :
701 : template <typename OutputType>
702 : void
703 316780 : MooseVariableFV<OutputType>::clearCaches()
704 : {
705 316780 : _elem_to_grad.clear();
706 316780 : }
707 :
708 : template <typename OutputType>
709 : unsigned int
710 5808 : MooseVariableFV<OutputType>::oldestSolutionStateRequested() const
711 : {
712 5808 : unsigned int state = 0;
713 5808 : state = std::max(state, _element_data->oldestSolutionStateRequested());
714 5808 : state = std::max(state, _neighbor_data->oldestSolutionStateRequested());
715 5808 : return state;
716 : }
717 :
718 : template <typename OutputType>
719 : void
720 31683 : MooseVariableFV<OutputType>::clearAllDofIndices()
721 : {
722 31683 : _element_data->clearDofIndices();
723 31683 : _neighbor_data->clearDofIndices();
724 31683 : }
725 :
726 : template <typename OutputType>
727 : typename MooseVariableFV<OutputType>::ValueType
728 86240534 : MooseVariableFV<OutputType>::evaluate(const FaceArg & face, const StateArg & state) const
729 : {
730 86240534 : const FaceInfo * const fi = face.fi;
731 : mooseAssert(fi, "The face information must be non-null");
732 86240534 : if (isDirichletBoundaryFace(*fi, face.face_side, state))
733 3306583 : return getDirichletBoundaryFaceValue(*fi, face.face_side, state);
734 82933951 : else if (isExtrapolatedBoundaryFace(*fi, face.face_side, state))
735 : {
736 678228 : bool two_term_boundary_expansion = _two_term_boundary_expansion;
737 678228 : if (face.limiter_type == Moose::FV::LimiterType::Upwind)
738 219 : if ((face.elem_is_upwind && face.face_side == fi->elemPtr()) ||
739 0 : (!face.elem_is_upwind && face.face_side == fi->neighborPtr()))
740 219 : two_term_boundary_expansion = false;
741 678228 : return getExtrapolatedBoundaryFaceValue(
742 678228 : *fi, two_term_boundary_expansion, face.correct_skewness, face.face_side, state);
743 : }
744 : else
745 : {
746 : mooseAssert(this->isInternalFace(*fi),
747 : "We must be either Dirichlet, extrapolated, or internal");
748 82255723 : return Moose::FV::interpolate(*this, face, state);
749 : }
750 : }
751 :
752 : template <typename OutputType>
753 : typename MooseVariableFV<OutputType>::ValueType
754 0 : MooseVariableFV<OutputType>::evaluate(const NodeArg & node_arg, const StateArg & state) const
755 : {
756 0 : const auto & node_to_elem_map = this->_mesh.nodeToElemMap();
757 0 : const auto & elem_ids = libmesh_map_find(node_to_elem_map, node_arg.node->id());
758 0 : ValueType sum = 0;
759 0 : Real total_weight = 0;
760 : mooseAssert(elem_ids.size(), "There should always be at least one element connected to a node");
761 0 : for (const auto elem_id : elem_ids)
762 : {
763 0 : const Elem * const elem = this->_mesh.queryElemPtr(elem_id);
764 : mooseAssert(elem, "We should have this element available");
765 0 : if (!this->hasBlocks(elem->subdomain_id()))
766 0 : continue;
767 0 : const ElemPointArg elem_point{
768 0 : elem, *node_arg.node, _face_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage};
769 0 : const auto weight = 1 / (*node_arg.node - elem->vertex_average()).norm();
770 0 : sum += weight * (*this)(elem_point, state);
771 0 : total_weight += weight;
772 : }
773 0 : return sum / total_weight;
774 0 : }
775 :
776 : template <typename OutputType>
777 : typename MooseVariableFV<OutputType>::DotType
778 : MooseVariableFV<OutputType>::evaluateDot(const ElemArg &, const StateArg &) const
779 : {
780 : mooseError("evaluateDot not implemented for this class of finite volume variables");
781 : }
782 :
783 : template <>
784 : ADReal
785 55560 : MooseVariableFV<Real>::evaluateDot(const ElemArg & elem_arg, const StateArg & state) const
786 : {
787 55560 : const Elem * const elem = elem_arg.elem;
788 : mooseAssert(state.state == 0,
789 : "We dot not currently support any time derivative evaluations other than for the "
790 : "current time-step");
791 : mooseAssert(_time_integrator && _time_integrator->dt(),
792 : "A time derivative is being requested but we do not have a time integrator so we'll "
793 : "have no idea how to compute it");
794 :
795 55560 : Moose::initDofIndices(const_cast<MooseVariableFV<Real> &>(*this), *elem);
796 :
797 : mooseAssert(
798 : this->_dof_indices.size() == 1,
799 : "There should only be one dof-index for a constant monomial variable on any given element");
800 :
801 55560 : const dof_id_type dof_index = this->_dof_indices[0];
802 :
803 55560 : if (_var_kind == Moose::VAR_SOLVER)
804 : {
805 55560 : ADReal dot = (*_solution)(dof_index);
806 104136 : if (ADReal::do_derivatives && state.state == 0 &&
807 48576 : _sys.number() == _subproblem.currentNlSysNum())
808 48576 : Moose::derivInsert(dot.derivatives(), dof_index, 1.);
809 55560 : _time_integrator->computeADTimeDerivatives(dot, dof_index, _ad_real_dummy);
810 55560 : return dot;
811 55560 : }
812 : else
813 0 : return (*_sys.solutionUDot())(dof_index);
814 : }
815 :
816 : template <>
817 : ADReal
818 0 : MooseVariableFV<Real>::evaluateDot(const FaceArg & face, const StateArg & state) const
819 : {
820 0 : const FaceInfo * const fi = face.fi;
821 : mooseAssert(fi, "The face information must be non-null");
822 0 : if (isDirichletBoundaryFace(*fi, face.face_side, state))
823 0 : return ADReal(0.0); // No time derivative if boundary value is set
824 0 : else if (isExtrapolatedBoundaryFace(*fi, face.face_side, state))
825 : {
826 : mooseAssert(face.face_side && this->hasBlocks(face.face_side->subdomain_id()),
827 : "If we are an extrapolated boundary face, then our FunctorBase::checkFace method "
828 : "should have assigned a non-null element that we are defined on");
829 0 : const auto elem_arg = ElemArg({face.face_side, face.correct_skewness});
830 : // For extrapolated boundary faces, note that we take the value of the time derivative at the
831 : // cell in contact with the face
832 0 : return evaluateDot(elem_arg, state);
833 : }
834 : else
835 : {
836 : mooseAssert(this->isInternalFace(*fi),
837 : "We must be either Dirichlet, extrapolated, or internal");
838 0 : return Moose::FV::interpolate<ADReal, FunctorEvaluationKind::Dot>(*this, face, state);
839 : }
840 : }
841 :
842 : template <>
843 : ADReal
844 0 : MooseVariableFV<Real>::evaluateDot(const ElemQpArg & elem_qp, const StateArg & state) const
845 : {
846 0 : return evaluateDot(ElemArg({elem_qp.elem, /*correct_skewness*/ false}), state);
847 : }
848 :
849 : template <typename OutputType>
850 : void
851 108528 : MooseVariableFV<OutputType>::prepareAux()
852 : {
853 108528 : _element_data->prepareAux();
854 108528 : _neighbor_data->prepareAux();
855 108528 : }
856 :
857 : template <typename OutputType>
858 : void
859 27035 : MooseVariableFV<OutputType>::determineBoundaryToDirichletBCMap()
860 : {
861 : mooseAssert(!Threads::in_threads,
862 : "This routine has not been implemented for threads. Please query this routine before "
863 : "a threaded region or contact a MOOSE developer to discuss.");
864 :
865 27035 : _boundary_id_to_dirichlet_bc.clear();
866 27035 : std::vector<FVDirichletBCBase *> bcs;
867 :
868 : // I believe because query() returns by value but condition returns by reference that binding to a
869 : // const lvalue reference results in the query() getting destructed and us holding onto a dangling
870 : // reference. I think that condition returned by value we would be able to bind to a const lvalue
871 : // reference here. But as it is we'll bind to a regular lvalue
872 27035 : const auto base_query = this->_subproblem.getMooseApp()
873 27035 : .theWarehouse()
874 : .query()
875 27035 : .template condition<AttribSystem>("FVDirichletBC")
876 27035 : .template condition<AttribThread>(_tid)
877 27035 : .template condition<AttribVar>(_var_num)
878 27035 : .template condition<AttribSysNum>(this->_sys.number());
879 :
880 120085 : for (const auto bnd_id : this->_mesh.getBoundaryIDs())
881 : {
882 93050 : auto base_query_copy = base_query;
883 186100 : base_query_copy.template condition<AttribBoundaries>(std::set<BoundaryID>({bnd_id}))
884 93050 : .queryInto(bcs);
885 : mooseAssert(bcs.size() <= 1, "cannot have multiple dirichlet BCs on the same boundary");
886 93050 : if (!bcs.empty())
887 38815 : _boundary_id_to_dirichlet_bc.emplace(bnd_id, bcs[0]);
888 : }
889 :
890 27035 : _dirichlet_map_setup = true;
891 27035 : }
892 :
893 : template <typename OutputType>
894 : void
895 27035 : MooseVariableFV<OutputType>::determineBoundaryToFluxBCMap()
896 : {
897 : mooseAssert(!Threads::in_threads,
898 : "This routine has not been implemented for threads. Please query this routine before "
899 : "a threaded region or contact a MOOSE developer to discuss.");
900 :
901 27035 : _boundary_id_to_flux_bc.clear();
902 27035 : std::vector<const FVFluxBC *> bcs;
903 :
904 : // I believe because query() returns by value but condition returns by reference that binding to a
905 : // const lvalue reference results in the query() getting destructed and us holding onto a dangling
906 : // reference. I think that condition returned by value we would be able to bind to a const lvalue
907 : // reference here. But as it is we'll bind to a regular lvalue
908 27035 : const auto base_query = this->_subproblem.getMooseApp()
909 27035 : .theWarehouse()
910 : .query()
911 27035 : .template condition<AttribSystem>("FVFluxBC")
912 27035 : .template condition<AttribThread>(_tid)
913 27035 : .template condition<AttribVar>(_var_num)
914 27035 : .template condition<AttribSysNum>(this->_sys.number());
915 :
916 120085 : for (const auto bnd_id : this->_mesh.getBoundaryIDs())
917 : {
918 93050 : auto base_query_copy = base_query;
919 186100 : base_query_copy.template condition<AttribBoundaries>(std::set<BoundaryID>({bnd_id}))
920 93050 : .queryInto(bcs);
921 93050 : if (!bcs.empty())
922 26773 : _boundary_id_to_flux_bc.emplace(bnd_id, bcs);
923 : }
924 :
925 27035 : _flux_map_setup = true;
926 27035 : }
927 :
928 : template class MooseVariableFV<Real>;
929 : // TODO: implement vector fv variable support. This will require some template
930 : // specializations for various member functions in this and the FV variable
931 : // classes. And then you will need to uncomment out the line below:
932 : // template class MooseVariableFV<RealVectorValue>;
|
